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Recent  work  has  demonstrated  the  power  of  combining  group  theory  with 
metaheuristic  search  methodologies  to  solve  discrete  optimization  problems.  Group 
theory  provides  tools  to  characterize  the  underlying  structures  in  move 
neighborhoods,  solution  representations  and  solution  landscapes.  Exploiting  these 
structures  with  group  theoretic  techniques  produces  highly  effective  and  efficient 
search  algorithms. 

Using  group  theory  we  develop  a  methodology  for  partitioning  the  solution 
space  into  orbits.  The  partitioning  is  performed  by  clustering  the  variables  based  on 
the  reduced  costs  of  the  LP  relaxation  creating  “good”  and  “bad”  orbits.  We  are  able 
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to  calculate  the  size  of  each  orbit  and  place  upper  and  lower  bounds  on  the  solutions 
contained  within.  The  search  efforts  can  then  be  directed  on  the  “good”  orbits. 

Based  on  these  ideas,  we  develop  a  Group  Theoretic  Tabu  Search  (GTTS) 
algorithm  for  solving  the  unicost  Set  Covering  Problem  (SCP),  GTTS-USCP.  We 
tested  our  algorithm  on  65  benchmark  problems  and  compared  the  results  against  the 
previous  best  known  and  solutions  obtained  by  CPLEX  9.0.  GTTS-USCP  discovered 
46  new  best  known  solutions.  GTTS-USCP  converged  significantly  faster  than 
CPLEX  for  all  problem  sets. 

We  explore  the  general  integer  linear  program  (ILP)  by  way  to  the  group 
minimization  problem  (GMP).  By  examining  the  local  search  in  terms  of  the  GMP, 
we  gain  insights  that  will  help  us  solve  the  ILP.  We  describe  the  local  search  for  the 
comer  polyhedron  in  the  space  of  the  non-basic  variables.  Integer  points  in  the  corner 
polyhedron  that  produce  an  all  integer  basis  form  a  sub-lattice.  We  develop  identity 
move  neighborhoods  that  allow  the  local  search  to  traverse  this  sub-lattice.  We  also 
develop  bound  strengthening  of  the  non-basic  variables  based  on  reduced  costs. 
These  bounds  effectively  shrink  the  corner  polyhedra  reducing  the  size  of  the  solution 
space  we  must  search. 

Based  on  this  research,  we  develop  a  GTTS  algorithm  for  solving  the  GMP, 
GTTS-GMP.  Since  the  GMP  can  be  formed  from  any  ILP,  this  algorithm  solves  the 
general  ILP.  The  algorithm  performs  well  on  a  diverse  set  of  benchmark  problems. 
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Chapter  1  -  Introduction 


Recent  work  has  demonstrated  the  power  of  combining  group  theory  with 
metaheuristic  search  methodologies  to  solve  discrete  optimization  problems.  Group 
theory  provides  tools  to  characterize  the  underlying  structures  in  move 
neighborhoods,  solution  representations  and  solution  landscapes.  Exploiting  these 
structures  with  group  theoretic  techniques  produces  highly  effective  and  efficient 
search  algorithms. 

Discrete  optimization  problems  may  be  divided  into  three  distinct  groups: 
partitioning,  ordering  and  partitioning-and-ordering  problems.  Partitioning  problems 
such  as  set  covering,  knapsack  and  min-cut  network  flow  problems  have  no  ordering 
context  and  require  only  that  the  solution  variables  be  placed  into  mutually  exclusive 
sets.  Ordering  problems  such  as  single-agent  traveling  salesman,  single-machine  job 
shop  scheduling  and  single-vehicle  routing  problems  require  that  a  permutation  of  the 
solution  variables  be  stipulated.  Partitioning-and-ordering  problems  such  as  multiple- 
agent  traveling  salesmen,  multiple-machine  job  shop  scheduling  and  multiple-vehicle 
routing  problems  require  that  the  solution  variables  be  partitioned  and  ordered  within 
each  partition. 

Since  previous  group  theoretic  metaheuristic  research  has  focused  on  either 
ordering  or  partitioning-and-ordering  problems,  only  the  symmetric  group  on  n 
letters,  S,„  the  group  of  permutations  of  n  distinct  objects,  has  been  employed.  S„  is 
inappropriate  for  strict  partitioning  problems.  An  appropriate  direction  for  current 
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research  is  to  investigate  whether  group  theory  can  improve  the  performance  of 
metaheuristic  search  methodologies  for  partitioning  problems. 

In  this  dissertation,  the  use  of  group  theory  for  the  characterization  of  move 
neighborhoods,  solution  representation  and  landscape  structures  for  partitioning 
problems  is  discussed.  Search  methodologies  that  exploit  these  structures  appear  to 
achieve  better  solutions  in  less  time  than  competing  approaches. 

In  addition,  Gomory’s  long  neglected  group  minimization  problem  (GMP)  is 
addressed.  The  GMP  seeks  to  construct  the  optimal  solution  to  an  integer  linear 
problem  (ILP)  from  the  optimal  solution  to  its  linear  relaxation.  In  the  1 970s,  there 
were  several  attempts  to  solve  the  GMP  with  dynamic  programming.  This  research 
reexamines  the  GMP  from  a  metaheuristic  perspective  and  describes  local  search 
neighborhoods  for  solving  the  general  ILP  in  this  context. 

We  will  now  overview  the  remaining  chapters  in  this  dissertation.  Chapter  2 
provides  a  brief  overview  of  group  theory  and  a  review  of  previous  operations 
research  literature  involving  group  theory.  Chapter  3  considers  the  general 
partitioning  problem  and  presents  group  theoretic  strategies  for  partitioning  the 
solution  space  to  enhance  the  exploration  of  that  space.  Chapter  4  provides  a  highly 
effective  and  efficient  reactive  tabu  search  (RTS)  implementation  of  one  of  these 
strategies  for  the  unicost  set  covering  problem  (SCP). 

Chapter  5  describes  move  neighborhoods  for  the  general  IP  in  terms  of  group 
theory  and  the  GMP.  Since  any  ILP  can  be  transformed  into  a  GMP,  these  ideas  can 
be  applied  withifi  most  metaheuristics  to  create  a  search  algorithm  that  solves  the 
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general  ILP.  Chapter  6  presents  a  powerful  RTS  algorithm  for  the  GMP  employing 
these  neighborhoods  and  ideas. 

Finally,  Chapter  7  provides  conclusions  from  this  research  as  well  as 
directions  for  further  research. 
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Chapter  2  -  Literature  Review 


This  chapter  presents  a  synopsis  of  group  theory  and  the  previous  research 
relevant  to  the  goals  of  this  dissertation.  Section  2.1  presents  a  brief  overview  of 
group  theory.  A  more  comprehensive  treatment  can  be  found  in  countless  abstract 
algebra  texts  such  as  Fraleigh  (1976)  or  Herstein  (1975).  Colletti  (1999)  provides  a 
robust  treatment  of  group  theory  from  the  perspective  of  metaheuristics.  Section  2.2 
provides  a  short  discussion  on  metaheuristics,  focusing  on  tabu  search  (TS).  Section 
2.3  provides  a  literature  review  on  the  use  of  group  theory  in  operations  research 
(OR).  Section  2.4  reviews  recent  work  on  group  theoretic  metaheuristic  search 
methods. 

2.1  An  Introduction  to  Group  Theory 


2.1.1  Groups 

Given  a  set  of  elements  G  and  a  binary  “multiplication”  operation  ©  then 
(G,©)  is  a  group  if 

1.  Vg,  h  gG  =>  g  ©  /z  e  G  (Closure) 

2.  g  ©  (h  ©  j)  =  (g  ©  h)  ©  j  \/g,h,j  e  G  (Associativity) 

3.3 e  eG  s.t.  g  © e  =  g  Vg  e  G  (Identity) 

4.  3g_1  e  G  s.t.  g  ©  g-1  =e  Vg  e  G  (Inverse) 

An  abelian  group  embodies  a  commutative  binary  operation  (Fraleigh  1976,  Herstein 
1975). 
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2.1.2  Subgroups 

Let  G  be  a  group  with  H  cG.  If  H  is  also  a  group  under  the  operation  ©  of 
G,  H  is  a  subgroup  of  G  denoted  H  <G .  H  must  be  closed  under  © ,  contain  the 
identity  element  e,  and  contain  g_1  eH  Vg  e  H  . 

2.1.3  Cosets 

If  H  <  Gandg  eG,  gH  ~{g®h,  he  H)  is  a  left  coset  of  H  in  G  and  Hg  is  a 
right  coset  of  H  in  G.  The  left  (right)  cosets  are  disjoint  and  partition  G  into  equal 
sized  sets.  The  number  of  left  (right)  cosets  is  |G|/|//|  (Colletti  1999).  If  G  is 
abelian,  gH  =  Hg  Vg  e  G . 

2.1.4  Cyclic  Groups 

For  n  >  0,  define  g"  as  the  n-fold  product  of  g  (i.e.,  g  multiplied  by  itself  n 
times).  For  n  <  0,  g"  is  the  n-fold  product  of  g'1 .  Finally,  define  g°  =  e.  The  order  of 
g  is  the  minimum  positive  value  of  i  such  that  g'  =  g°  =  e. 

Group  closure  implies  if  g  eH  then  g'  eH  \He7L.  Let  g  eG  then  the  cyclic 

group,  H  ={gn\ne%  }  =  (g)  is  the  smallest  subgroup  of  G  that  contains  g 

(Fraleigh  1976).  The  group  K  =  (g,h,j)  is  the  smallest  subgroup  that  contains  g,  h, 
and  j. 
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2.1.5  The  External  Direct  Product  of  Groups 

Let  G\,  G2,  Gn  be  a  collection  of  groups.  Define  the  group  G  to  be  the 
Cartesian  product  of  the  G„  i.e.,  the  external  direct  product  of  the  groups  G;.  The 
binary  operation  in  G  is  component-wise  multiplication: 

{gvg2’---’8n){KK---A)  =  {gA’gi}h’---’gnK)-  The  identity  of  G  is 
(e1,e2,...,en)and  the  inverse  of  any  element  is  (g,,g2,...,gB)-1 
(Fraleigh  1976,  Herstein  1975).  If  all  G \  are  abelian  then  G  is  abelian. 

2.1.6  Factor  Groups 

For  g,heG,  define  gh  =  h~lgh  to  be  the  conjugate  of  g  by  h  (Fraleigh 
1976).  If  3 jeG  such  that  gJ=h,  g  and  h  are  conjugates.  For  abelian  G, 
gh  =  h~'gh  =  hh~'g  =  g  and  Vg  e  G  the  only  conjugate  element  is  g  itself.  N  <  G  is 
a  normal  subgroup  of  G,  denoted  N  <G ,  if  Ng  =  N  Vg  e  G  .  For  abelian  G,  all 
subgroups  are  normal. 

If  N  <  G ,  the  left  cosets,  gN  Vg  e  G ,  form  the  factor  group  of  G  modulo  N , 
denoted  GIN,  under  the  set  product  operation  (Fraleigh  1976): 

A®  B  ={ab :  a  e  A,b  e  B]VA,B  c  G 
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2.1.7  Homomorphism 

A  function  /  mapping  between  two  possibly  distinct  groups,  (G,,©)  and 
(G2,®) ,  is  a  homomorphism  if  the  mapping  preserves  the  group  operation.  If 
f:Gl-^G2  then /is  a  homomorphism  iff  f  (g  ©  h)  =  f  (g)®  f  (h)  Vg,/jeG,.  If  in 

addition  the  homomorphism  is  a  bijection  then  it  is  an  isomorphism.  Two  groups 
which  are  isomorphic  are  essentially  identical  with  the  elements  relabeled.  If  N  <  G 
then  the  natural  homomorphism  is  f  :  G  — >•  G  /  N  mapping  each  g  eG  into  the 
coset  gN  (Colletd  1999). 

2.1.8  Group  Action 

Given  a  group  G  and  a  set  T,  the  group  action  of  G  on  T,  gT,  is  a  remapping 
of  T  such  that  Vg  e  G  and  Vs  e  T ,  we  have  sAg=teT  (Colletti  1999).  Additionally 
gT  must  have  the  following  properties: 

1 .  tAe  =  t  Vt  eT  (where  e  is  the  identity  of  G) 

2.  (tAg)Ah  =  t A(gh)  VteT  and g,heG 

A  group  action  partitions  T  into  disjoint  orbits.  For  group  action  qT,  the  orbit  of 
teT  is  Orbit(G,t)  =  {tAg|g  eGj.  If  s,teT  are  in  the  same  orbit  then  3  geG 

such  that  sAg=t. 
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2.1.9  Examples 


Sn  has  been  used  in  significant  recent  work  in  metaheuristics  (Colletti  1999, 
Crino  2002,  Wiley  2001).  A  permutation  may  be  represented  by  a  2  by  n  matrix 
corresponding  to  a  1-1  and  onto  mapping  of  the  integers  { 1,2, ...,«}  where  the  integer 
in  the  top  row  is  replaced  in  the  order  by  its  image  in  the  bottom  row.  For  example,  if 
n=6  then 


^1  2  3  4  5  6" 

,3  1  2  4  6  5, 


"1  2  3  4  5  6" 

,6  1  5  3  2  4, 


p  and  q  are  permutations.  The  set  of  all  permutations  of  degree  n  along  with  the 
binary  operation  of  function  composition,  p®q=pq  =  q(p(x)) ,  form  Sn . 

Permutations  are  often  written  in  cycle  notation  for  example  p  =  (132)(4)(56) 
=  (132)(56).  Single  character  cycles  or  1 -cycles  map  a  letter  onto  itself  and  are 
typically  not  shown.  Cycles  with  2  characters  (2-cycles)  are  transpositions.  The 
letters  moved  by  p  are  denoted  move(p).  In  chapters  3  and  4,  we  will  use  a  group 
action  based  on  a  direct  product  of  Sn . 


Another  important  group  in  this  research,  G(D ) ,  is  based  on  fractions  under 
addition  modulo  1.  The  group  elements  are  k/D  for  some  k,D  e 7L,  0 <k<D. 
G(D )  is  closed  under  addition  modulo  1,  has  identity  0/D,  and  each  element  k/D  has 


inverse  ( D-k )  /  D.  G(D)  is  abelian,  cyclic  with  generator  1/D  and  has  size  D.  Since 

we  can  relabel  the  elements  by  multiplying  each  by  D,  the  group  is  isomorphic  to  the 
set  of  integers  under  addition  modulo  D,  denoted  Z(D).  The  group  Z(D)  consists  of 
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the  integers  {0,  1...  D-l}.  In  chapters  5  and  6,  we  will  use  a  subgroup  of  an  external 
direct  product  of  m  of  these  groups. 

2.2  HeuristidVIethods 

Many  discrete  optimization  problems  are  NP-hard  (Wolsey  1998).  Such 
problems  are  difficult,  i.e.,  polynomial  time  algorithms  to  solve  these  problems  do  not 
exist.  For  such  problems,  the  number  of  solutions  grows  exponentially  with  problem 
size  and  all  solutions  must  be  explicitly  or  implicitly  enumerated  to  guarantee 
optimality.  Consequently,  heuristic  methods  are  often  used  to  address  these  types  of 
problems.  While  no  guarantee  of  solution  quality  is  achieved,  empirical  results  have 
shown  these  methods  often  return  acceptable  solutions  in  a  very  short  amount  of  time. 

2.2.1  Local  Search  Methods 

Heuristic  methods  such  as  steepest  descent,  simulated  annealing  and  tabu 
search  are  local  search  methods  (Glover  and  Laguna  1997,  Reeves  1995).  In  a  local 
search  method,  a  starting  solution  is  chosen.  A  move  modifies  the  current  incumbent 
solution  and  all  solutions  that  can  be  reached  in  one  move  comprise  the  incumbent’s 
neighborhood.  The  new  solution  is  chosen  from  the  neighborhood  based  on  a  merit 
function  and  the  details  of  the  algorithm.  The  new  solution  becomes  the  incumbent 
and  the  procedure  repeats  until  a  termination  criterion  is  satisfied. 
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2.2.2  Tabu  Search 


In  tabu  search  (TS)  (Glover  and  Laguna  1997),  we  prohibit  recently  visited 
solutions  from  being  revisited  for  tabu  tenure  iterations.  At  each  iteration  the  best 
non-tabu  neighboring  solution  is  selected.  The  tabu  memory  structure  allows  escape 
from  local  optima  to  continue  the  search.  TS  has  been  shown  to  be  quite  effective  in 
solving  complex  optimization  problems  (Crino  2002,  Glover  and  Laguna  1997,  Wiley 
2001). 

Reactive  Tabu  Search  (RTS)  was  developed  by  Battiti  and  Tecchiolli  (1994). 
In  early  TS  implementations,  the  length  of  the  tabu  tenure  was  fixed  or  randomly 
chosen.  In  RTS,  the  current  tabu  tenure  is  based  on  the  occurrence  of  repeated  visits 
to  the  same  solutions.  When  the  algorithm  detects  sufficient  repetitions,  the  tabu 
tenure  is  increased  to  encourage  the  algorithm  to  diversify  into  a  region  of  the 
solution  space  not  yet  explored. 

The  tabu  memory  structure  in  RTS  also  helps  detect  the  occurrence  of  cycling 
-  revisiting  the  same  solutions  repeatedly.  When  the  algorithm  detects  cycling,  an 
escape  procedure  is  implemented  to  attempt  to  break  the  cycle  and  hopefully  diversify 
into  a  region  of  the  solution  space  not  yet  explored. 

2.2.3  Landscape  Theory 

An  often  ignored  but  very  important  aspect  of  local  search  methods  is  the 
problem  landscape.  The  landscape  structure  is  determined  by  the  objective  function, 
the  solution  space  and  the  neighborhood  definition  (Barnes  et  al.  2003).  The 
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neighborhood  definition  determines  whether  the  solution  landscape  will  be  favorable 
for  local  search. 

Weinberger  (1990)  addresses  the  idea  of  correlated  and  uncorrelated 
landscapes.  He  classifies  landscapes  with  approximate  decaying  exponential 
autocorrelation  spectra  to  be  AR(1)  (first  order  autoregressive)  landscapes. 
Considering  only  regular  and  symmetric  neighborhoods,  Grover  (1992)  derives  a 
difference  equation  for  the  neighborhoods  similar  to  a  well  known  differential 
equation  used  in  mathematical  physics.  For  landscapes  that  satisfy  Grover’s 
difference  equation,  all  local  minima  are  less  than  the  average  solution  value  over  the 
landscape  and  all  local  maxima  are  greater  than  the  average  solution  value. 

Using  group  theory,  Colletti  (1999)  proves  that  symmetric  multiple  traveling 
salesmen  problems  (m-STSPs)  with  move  neighborhoods  based  on  an  arbitrary 
collection  of  k-city  exchange  moves  will  satisfy  Grover’s  difference  equation  for  any 
value  of  k. 

Again  considering  only  regular  and  symmetric  neighborhoods  Stadler  (1996) 
develops  his  Laplacian  and  forms  the  matrix  version  Grover’s  equation.  He  shows 
that  if  the  normalized  objective  vector  is  an  eigenvector  of  his  Laplacian  matrix  then 
the  landscape  will  satisfy  Grover’s  equation.  He  calls  the  associated  landscapes 
elementary.  He  also  claims  that  for  regular  symmetric  transition  matrices,  the 
associated  landscape  is  elementary  if  and  only  if  it  is  an  AR(1)  landscape. 

Barnes  et  al.  (2003)  extend  the  theory  of  elementary  landscapes  by  using  a 
more  general  Laplacian  matrix  that  also  embraces  arbitrary  transition  matrices.  They 
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also  define  two  types  of  elementary  landscapes:  smooth  and  rugged.  Dimova  et  al. 
(2005)  extend  Stadler’s  results  on  AR(1)  processes  to  arbitrary  transition  matrices. 
For  any  transition  matrix,  the  landscape  is  elementary  if  and  only  if  the  random  walk 
on  the  landscape  follows  an  AR(1)  process. 

2.3  Group  Theory  and  Operations  Research 

2.3.1  Integer  Programming  and  Cutting  Planes 

Gomory  (1963,  1965,  1967,  1969)  developed  a  methodology  for  adding  valid 
inequalities  or  cutting  planes  to  an  integer  linear  program  (ILP)  to  drive  the  relaxed 
linear  program  (LP)  to  an  integer  solution.  The  inequalities  are  formed  using  the 
fractional  components  of  an  integer  combination  of  the  constraint  rows  in  the  LP 
optimal  tableau.  Each  inequality  renders  the  current  LP  optimum  infeasible.  In 
principle,  repeatedly  adding  cuts  and  resolving  will  eventually  yield  the  optimal 
solution  to  the  original  ILP. 

Gomory  (1963)  shows: 

(1)  The  fractional  components  take  the  form  k/D  for  some  k  e  ,  0 <k<D 
where  D  is  the  determinant  of  the  optimal  basis  to  the  original  LP 
relaxation. 

(2)  The  fractional  cuts  form  a  group  under  component-wise  addition  mod  1 . 
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(3)  The  group  is  isomorphic  to  a  factor  group  M(/)/M(B)  where  M(/)  is  the 
group  of  all  integer  vectors  of  size  m  and  M(B)  is  the  subgroup  of  all 
linear  integer  combinations  of  the  columns  in  the  optimal  LP  basis. 

Pure  cutting  plane  methodologies  fell  out  of  favor  in  the  1980s  due  to  their 
lack  of  convergence  in  practical  implementations.  Nevertheless,  the  combination  of 
cutting  planes  and  branch  and  bound  algorithms  (branch  &  cut)  are  among  the  most 
popular  and  effective  methodologies  for  solving  ILPs  and  mixed  ILPs  today. 

2,3.2  The  Group  Minimization  Problem 

The  group  minimization  problem  (GMP)  (Gomory  1965,  1967,  1969)  is  a 
mathematical  formulation  for  solving  ILPs  (Johnson  1980)  by  perturbing  the  non- 
basic  variables  from  the  associated  LP  relaxation  while  ensuring  the  integrality  of  the 
basic  variables.  Since  any  such  perturbation  degrades  solution  quality,  the  goal  is  to 
minimize  such  changes.  When  the  GMP  yields  non-negative  basic  variables,  it  solves 
the  original  ILP,  i.e.,  an  algorithm  that  solves  the  GMP  also  solves  the  general  ILP. 

The  GMP  can  be  expressed  using  either  one  of  two  abelian  groups,  H{a )  or 
M(I)/M(B).  As  alluded  to  in  Section  2.1.9,  H(a)  is  a  subgroup  of  an  external  direct 
product  of  m  G(D )  groups  where  m  is  the  number  of  rows  in  the  basis.  The  group 

operation  is  addition  modulo  1  and  the  group  is  generated  by  the  fractional 
components  of  the  non-basic  columns  in  the  optimal  LP  tableau.  The  size  of  the  full 
group  containing  H{a)  is  Dm. 
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M(/)/M(B)  is  the  factor  group  described  above.  Each  non-basic  column  is 
mapped  to  an  element  of  the  factor  group  using  the  natural  homomorphism.  The 
elements  then  generate  the  full  factor  group  which  also  is  of  size  <  D.  M(7)/M(2J)  is 
isomorphic  to  H(a). 

Gomory  (1969)  shows  how  the  GMP  can  be  used  to  transform  the  problem 
into  an  integer  program  with  one  constraint.  The  problem’s  decision  variables  are  the 
integer  multipliers  for  each  group  element.  Early  GMP  solution  approaches  (Glover 
1968,  Shapiro  1968a,  1968b)  attempted  to  enumerate  the  group  using  dynamic 
programming  or  network  flow  optimization.  Since  these  approaches  could  manage 
only  trivially  sized  GMPs,  the  GMP  was  largely  forgotten  and  subsequently  ignored 
by  the  operations  research  community. 

2.3.3  Corner  Polyhedra 

Gomory  (1967,  1969)  described  the  geometry  associated  with  his  cutting 
planes  and  GMP.  Let  P  be  the  feasible  region  for  the  LP  relaxation  of  an  ILP.  Given 
any  vertex  of  P,  we  relax  all  constraints  not  passing  through  the  vertex  to  form  a 
cone,  K.  The  convex  hull  of  the  integer  points  in  K  is  the  corner  polyhedra  (CP)  for 
that  vertex. 

Gomory  shows  if  the  vertex  is  LP  optimal  then  the  optimal  solution  to  the  ILP 
is  a  vertex  on  the  corresponding  CP.  Gomory’s  fractional  cuts  are  the  faces  of  CP. 
CP  can  be  expressed  in  terms  of  the  original  decision  variables  or  in  terms  of  the  non- 
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basic  variables.  The  feasible  region  for  the  GMP  corresponds  to  the  CP  in  terms  of 
the  non-basic  variables. 

Corner  polyhedra  may  also  be  expressed  in  terms  of  the  full  group  generated 
by  the  elements  associated  with  the  non-basic  columns.  This  increases  the  dimension 
of  CP  to  D.  This  is  the  master  CP.  Any  facet  of  the  CP  is  a  facet  of  the  master  CP 
and  any  facet  of  the  master  CP  is  a  face  of  the  CP  but  not  necessarily  a  facet.  Many 
problems  with  different  CP  will  share  a  common  master  CP.  Gomory  (1969) 
generated  facets  of  the  master  CP  in  hopes  of  finding  facets  to  all  corresponding  CP. 

2.3.4  Recent  Work 

Gomory’ s  work  with  corner  polyhedra  and  master  corner  polyhedra  to 
generate  cutting  planes  has  mostly  lain  dormant  for  the  last  two  decades.  Recently, 
Gomory  and  others  have  begun  to  resurrect  the  research  (Gomory  and  Johnson  2003a, 
2003b,  Gomory  et  al.  2003,  Araoz  et  al.  2003).  The  focus  of  that  new  research  is  still 
on  cutting  planes  and  the  faces  of  the  CP.  The  research  documented  in  this 
dissertation  focuses  on  the  interior  of  the  CP. 

2.4  Group  Theory  and  Metaheuristics 

Although  not  commonly  acknowledged,  many  aspects  of  metaheuristics  can 
be  defined  in  terms  of  group  theory.  Colletti  (1999)  gives  a  comprehensive  treatment 
of  group  theory  in  the  context  of  metaheuristics.  Using  the  TSP  and  TS,  he  classifies 
and  defines  current  move  definitions  in  terms  of  group  theory  and  develops  composite 
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move  strategies  that  would  be  difficult  to  generate  using  other  methods.  He  provides 
efficient  methods  for  escaping  from  chaotic  attractors  during  the  search  and  for 
generating  search  neighborhoods. 

This  Group  Theoretic  Tabu  Search  (GTTS)  approach  has  been  very 
successfully  applied  to  several  complex  problems.  Wiley  (2001)  implements 
Colletti’s  ideas  to  solve  the  Aerial  Fleet  Refueling  Problem.  A  solution  representation 
is  developed  using  S„  and  dynamic  move  neighborhoods  are  developed  based  on 
conjugation  and  function  composition  to  explore  the  solution  space.  Crino  (2002) 
extended  these  ideas  in  solving  the  military  theater  distribution  problem.  He  used 
group  actions  to  partition  the  solution  space  into  orbits  and  orbital  planes.  Our  work 
in  chapters  3  and  4  is  based  on  Crino’s  approach. 

The  next  chapter  considers  the  general  partitioning  problem  and  presents  a 
group  theoretic  strategy  for  partitioning  the  solution  space  to  enhance  the  exploration 
of  that  space. 
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Chapter  3  -  Partitioning  the  Solution  Space 


Since  the  solution  space  for  a  typical  NP-hard  problem  is  immense, 
appropriate  partitioning  of  the  solution  space  can  improve  the  search  by  aiding  in 
diversification,  intensification,  and  cycle  prevention.  Group  Theory  provides  many 
methods  for  such  partitions. 

3.1  Partitioning  Into  Orbits  (Ordering  Problems) 

3.1.1  Theory 

Colletti  (1999)  represents  solutions  to  the  single  and  multiple  agent  TSPs  as 
elements  of  Sn  and  proposes  the  use  of  group  actions  and  orbits  to  partition  the 
solution  space.  For  example,  a  solution  to  a  10-city  2-TSP  may  be  (4  6  2  9)(3  1  10  8 
5  7)  with  agent  1  visiting  cities  4,  6,  2,  and  9  in  order  and  agent  2  visiting  cities  3,  1, 
10,  8,  5,  and  7.  Each  agent’s  tour  is  a  cycle  and  the  number  and  size  of  the  cycles  in 
the  solution  is  the  solution’s  cycle  structure. 

Two  elements  with  the  same  cycle  structure  are  conjugates  and  are  in  the 
same  conjugacy  class.  A  group  action  on  Sn  based  on  conjugation  with  H  <  S„  does 
not  change  cycle  structure.  A  conjugacy  class  with  b  a-cycles  and  d  c-cycles  is 
denoted  by  abcd. 
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Example  3.1 


Let  Sio  be  our  solution  space  and  let  H  =^(1  2),(l  3)^  ( H  <SI0).  Let  the  current 

solution  be  x  =  (4  6  2  9)(3  1  10  8  5  7)  £  4161  .  The  orbit  neighborhood  for  H  given  x, 
a'h  ,  is  all  solutions  reachable  by  rearranging  1,  2,  and  3  in  x,  i.e,  all  solutions 
constructed  by  conjugating  x  with  any  element  of  H: 

xA(l  2)  =  (4  6  1  9)(3  2  10  8  5  7) 
jca(1  3)  =  (4  6  2  9)(1  3  10  8  5  7) 
xA(2  3)  =  (4  6  3  9)(2  1  10  8  5  7) 
jcA(l  2  3)  =  (4  6  1  9)(2  3  10  8  5  7) 
jca(1  3  2)  =  (4  6  3  9)(1  2  10  8  5  7) 
xAe  =  (462  9)(3  1  10857) 

To  change  to  another  orbit  with  the  same  cycle  structure  we  can  use 
conjugation  with  a  permutation  not  in  the  subgroup  H.  We  must  use  some  other 
operation,  like  function  composition,  to  move  to  an  orbit  with  a  different  cycle 
structure.  A  transversal  of  the  orbits  is  a  list  containing  one  element  from  each 
possible  orbit.  Changing  the  current  solution  to  an  element  in  a  different  orbit  moves 
the  search  to  that  orbit  and  conjugation  with  the  subgroup  H  allows  the  search  to 
explore  the  orbit.  Unfortunately,  creating  transversals  can  be  quite  costly  as  the 
number  of  orbits  increases. 

Performing  the  group  action  on  Sn  using  all  of  Sn  partitions  the  solution  space 
into  orbits  each  representing  a  different  cycle  structure  or  conjugacy  class  (i.e.  1G1, 
1  !9x,  2181,  25,  3222,  etc.)  .  For  example,  orbit  456x  would  contain  all  solutions 
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assigning  4  cities  to  agent  1  and  6  cities  to  agent  2.  Alternatively,  we  can  perform  the 
group  action  using  subgroups  of  Sn. 

Colletti  (1999)  also  discusses  the  pros  and  cons  of  using  different  sized 
subgroups  of  Sn  for  the  group  action.  A  small  subgroup  yields  smaller  orbits  (each 
possibly  small  enough  to  search  exhaustively)  but  a  greater  number  of  such  orbits. 
Using  too  many  orbits  demands  an  effective  method  to  search  for  an  orbit  to  explore. 
A  larger  subgroup  yields  fewer  orbits,  likely  so  large  that  exploring  any  orbit  is 
difficult. 

3.1.2  Application 

Crino  (2002,  2004)  implements  Colletti ’s  ideas  to  solve  the  theater 
distribution  routing  and  scheduling  problem  (TDVRSP).  He  partitions  the  solution 
space  into  orbits  and  further  partitions  the  orbits  into  sub-orbits.  This  is  accomplished 
by  using  a  group  action  on  each  orbit  with  a  subgroup  of  the  group  used  to  create  the 
orbit.  Crino  uses  a  reactive  tabu  search  procedure  to  search  for  orbits  which  are  small 
enough  to  be  exhaustively  enumerated.  Once  an  orbit  has  been  explored  the  orbit  is 
made  tabu. 

The  first  partition  is  based  solely  on  cycle  structure  or  conjugacy  class.  The 
second  partition  is  based  on  cycle  order  (i.e.  2181  vs  8*2*  ).  Group  theory  does  not 
distinguish  elements  of  S„  by  cycle  order  abcd  is  equivalent  to  cdab;  however,  this 
distinction  is  required  as  the  TDVRSP  vehicles  assigned  to  each  cycle  are  non- 
homogeneous. 


19 


If  the  search  discovered  a  good  solution  with  4  cities  assigned  to  the  first 
vehicle  and  6  assigned  to  the  second,  it  is  likely  that  there  are  other  good  solutions  in 
the  same  orbit.  Similarly  if  the  search  has  trouble  finding  a  feasible  solution  in  orbit 
9' l1,  we  may  begin  to  believe  that  the  orbit  does  not  contain  any  feasible  solutions 
and  we  can  abandon  it.  While  we  may  not  know  which  orbits  are  good  or  bad  a 
priori,  we  may  be  able  to  narrow  down  the  good  orbits  shortly  into  the  search  and 
focus  our  efforts  there. 

Crino  further  partitions  the  above  orbits  by  assigning  customers  into  sets.  The 
number  and  size  of  the  sets  creates  the  first  partition,  an  orbital  plane.  The  orbital 
planes  are  then  partitioned  into  orbits  based  on  which  customers  are  assigned  to  each 
set.  All  of  these  partitions  are  created  by  a  group  action  with  some  subgroup  of  Sn. 
Exchanging  2  or  more  customers  between  the  sets  will  move  to  a  new  orbit  on  the 
orbital  plane.  Modifying  the  size  or  number  of  sets  will  move  to  a  new  orbital  plane. 

3.1.3  Observations 

Tabu  search  and  other  local  search  heuristics  have  been  shown  to  be  quite 
effective  at  finding  good  solutions  quickly  without  partitioning  (Glover  and  Laguna 
1997,  Combs  2004).  When  we  partition  the  solution  space,  we  restrict  the  movement 
of  the  search  algorithm.  If  this  is  a  good  restriction,  limiting  the  search  to  the  “good” 
areas  of  the  solution  space,  the  partitioning  can  be  quite  effective.  However,  arbitrary 
partitioning  or  too  many  partitions  may  trap  the  search  in  bad  regions  of  the  solution 
space. 
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The  first  two  levels  of  partitioning  used  by  Crino,  cycle  structure  and  cycle 
order,  have  some  relevance  to  the  problem  being  solved.  It  is  reasonable  to  assume 
that  this  partitioning  scheme  would  create  partitions  containing  solutions  of  similar 
value.  The  third  and  fourth  levels  of  partitioning  created  by  grouping  the  customers 
into  arbitrary  sets  have  no  logical  link  to  the  problems  being  solved.  One  good 
solution  in  an  orbital  plane  need  not  indicate  a  higher  or  lower  probability  of  other 
good  solutions  in  that  plane.  Superior  results  could  probably  be  achieved  if  only  the 
first  two  levels  of  partitioning  were  applied.  In  general,  partitioning  the  solution 
space  into  “good”  and  “bad”  partitions  will  allow  the  search  to  focus  on  the  good 
partitions. 

3.2  Partitioning  Into  Orbits  (Partitioning  Problems) 

For  strict  partitioning  problems,  using  Sn  for  the  solution  representation  is 
inappropriate.  A  vector  of  integers  suffices.  For  binary  programs,  the  solution 
representation  is  a  binary  vector  of  size  n  and  the  solution  space,  X,  is  the  set  of  all 
binary  vectors  of  size  n. 

We  can  create  a  group  action  on  X  using  Sn  or  a  subgroup  of  S„  and  partition 
the  solution  space  into  orbits.  Given  a  permutation  p  eSn  and  a  solution  vector  x,  we 

define  the  group  action  of  p  on  x  as  moving  the  value  at  vector  element  i  to  vector 
element  j  if  j  follows  i  in  p.  For  example,  if  *  =  (1012501310)  and  p  =  (  2  3  4  ) 
then  xAp  =  (1201501  3  1  0). 
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This  operation  is  a  valid  group  action  because  the  result  of  the  action  is  in  X 
and  the  properties  from  2.1.8  hold.  Using  all  of  Sn  for  the  group  action  on  X  creates 
only  a  single  orbit.  In  the  next  section,  we  show  how  to  use  a  subgroup  of  Sn  by 
clustering  the  elements  of  the  solution  and  restricting  movement  of  the  values  to 
within  the  clusters. 

3.2.1  Variable  Clustering 

Grouping  the  problem  variables  into  clusters  and  using  an  external  direct 
product  of  Sn  subgroups  creates  orbits  and  partitions  the  solution  space.  Let  there  be 
nk  variables  in  cluster  k  and  let  r  be  the  number  of  clusters.  The  group  acting  upon 

the  set  of  n  vectors  is  x...x Snr  where  ^=] nr  =  n .  For  example,  if  n  =  10,  n\  -  3, 

«2  =  3,  «3  =  4  and  r  =  3  a  solution  may  be  x  =  [l0l|250|l310].  The  group  acting 

upon  x  is  Sn3xSn3xSni  and  the  orbit  containing  x  contains  all  solutions  where  the 

first  cluster  contains  two  Is  and  a  0,  the  second  contains  a  2,  5,  and  0,  and  the  last 
contains  two  Is,  a  3  and  a  0.  The  orbit  of  x  is 

[l0l|250|l310]  [110|250|1310]  [01l|250|l310]  [l0l|520|l310]  [ll0|520|l310] 

[101|502|1310]  [110|502|1310]  [01 1 1502  |l  3 10]  [l0l|520|l310]  [01 1 1205  |l3 10] 

[lOl  |025  |l  3 10]  [l  1 0  [025  |l  3 1 0]  [01l|025|l310]  [l0l|205|l310]  [ll0|205|l310] 

[110|052|1130]  [ll0|052|3110]  [01l|052|3110]  [l0l|052|3110]  ... 
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If  we  define  cardinality  as  the  number  of  times  a  value  appears  in  the  cluster, 
an  orbit  can  be  characterized  by  the  cardinality  Of  the  components  of  each  cluster.  If 
a  problem  is  of  order  d  (i.e.  variable  values  are  0  thru  d- 1)  and  has  r  clusters,  we  can 
capture  the  orbit  characterization  in  a  d-1  by  r  matrix,  Op,  where  element  (ftk  is  the 
number  times  non-zero  value  i  appears  in  cluster  k  in  orbit  p.  For  example,  assume 
the  orbit  above  is  from  an  order  10  problem,  the  orbit  matrix  is 

f2  0  2l 
0  1  0 
0  0  1 
0  0  0 
op  =  0  1  0 
0  0  0 
0  0  0 
0  0  0 
[0  0  0] 

3.2.2  Which  Variables  Do  We  Cluster? 

Clustering  the  variables  is  easy.  Since  the  goal  is  to  group  solutions  of  similar 
quality  in  the  same  orbits,  determining  how  many  clusters  to  create  and  which 
variables  to  assign  to  each  cluster  can  be  more  challenging. 

Often,  a  logical  clustering  of  the  variables  is  obvious  from  the  context  of  the 
problem.  For  knapsack  problems,  it  seems  reasonable  to  assume  that  orbits  with 
more  non-zero  values  in  clusters  with  high  benefit/cost  ratios  would  be  superior  to 
other  orbits.  If  our  problem  contains  special  ordered  sets  (SOSs)  then  they  provide  a 
very  natural  clustering.  For  example,  if  we  have  3  SOSs  of  type  1  (only  1  member  of 
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the  set  can  be  non-zero),  we  create  3  clusters,  1  for  each  set.  Any  orbit  containing 
more  than  1  non-zero  element  in  each  cluster  is  an  infeasible  orbit  and  can  be 
discarded. 

Borrowing  from  the  ideas  of  Gomory  (1963,  1965,  1967,  1969),  we  can  solve 
the  LP  relaxation  of  our  problem  and  cluster  the  variables  based  on  reduced  costs. 
We  can  place  all  of  the  basic  variables  in  the  same  cluster  and  cluster  the  non-basic 
variables  separately  based  on  their  reduced  costs. 

We  cannot  use  reduced  costs  to  explicitly  determine  which  non-basic 
variables  will  change  values  in  the  near  optimal  solutions  (see  Section  3.3.3). 
However,  we  can  use  reduced  costs  as  a  heuristic  indicator  of  such  changes.  To 
reduce  the  number  of  orbits,  we  place  the  variables  in  the  same  cluster  if  their  reduced 
costs  are  within  a  certain  range  of  each  other. 

The  number  of  clusters  must  also  be  decided  and  may  also  be  determined  by 
the  problem  context.  If  it  is  not  clear  how  many  clusters  to  create,  it  is  preferable  to 
have  too  few  as  opposed  to  too  many.  Creating  too  many  clusters  leads  to  too  many 
orbits,  and  too  many  orbits  restricts  the  movement  of  the  search  making  it  difficult  to 
find  the  good  regions  of  the  solution  space.  If  we  assign  each  variable  to  its  own 
cluster,  exhaustive  enumeration  is  implied.  Alternatively,  if  we  assign  all  variables  to 
the  same  cluster,  general  tabu  search  is  implied.  Clearly  the  latter  is  preferable. 
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3.3  Orbit  Exploration 


3.3.1  Intra-Orbit  Move  Neighborhoods 

If  the  orbits  are  small  enough  we  may  be  able  to  enumerate  them;  however, 
this  is  not  likely  to  be  the  case.  So  we  must  create  a  move  neighborhood  that  is  of 
manageable  size,  restricts  the  search  to  the  current  orbit,  and  provides  connectivity  to 
all  solutions  within  the  orbit.  Let  G  =  Sn]x...xSnr  be  the  group  used  to  create  the 
orbits,  we  can  use  any  set  of  elements  H  cG  that  generates  G  as  a  move 
neighborhood  possessing  these  desired  properties.  H  is  a  subset  but  not  a  subgroup  of 
G.  A  subgroup  by  definition  is  closed  and  therefore  no  proper  subgroup  can  generate 
the  full  group,  therefore  H  cannot  be  a  subgroup  of  G. 

Every  element  of  H  is  in  G  and  since  G  is  closed  all  combination  of  elements 
of  H  are  also  in  G.  By  property  2  of  a  group  action,  we  know 
(jca  g) A  h  =  x*(gh)  Vxe  X  and  g,h  e  H  ,  so  executing  multiple  moves  in 

succession  is  equivalent  to  executing  one  move  using  the  product  of  the  move 
elements.  Therefore  the  neighborhood  based  on  H  restricts  the  search  to  the  current 
orbit.  Given  any  solution  in  the  orbit,  any  other  solution  in  the  orbit  can  be  generated 
using  an  element  of  G.  Since  the  elements  of  H  generate  G,  executing  multiple 
moves  will  eventually  generate  all  of  G  and  the  entire  orbit.  The  size  of  the 
neighborhood  is  simply  the  size  of  H. 
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Every  permutation  is  a  product  of  transpositions  and  so  the  set  of  all 


transpositions  in  S„k  generates  Snk.  For  G  the  number  of  transpositions  is 


This  set  is  essentially  a  within-cluster  swap  neighborhood  (including  duplicate 
solutions).  The  size  of  a  full  neighborhood  without  regard  to  the  clusters  would  be  a 

( 

much  larger  .  So  in  addition  to  partitioning  the  space  the  orbit  restriction  acts  as 

w 

a  candidate  list  restricting  neighborhood  size. 

There  are  even  smaller  neighborhoods  available.  The  subset  of  transpositions 
containing  all  transpositions  in  S„k  with  the  same  first  element  will  also  generate  S„k 
(i.e.  {(1  2),(l  3),...,(1  n,)})  as  will  the  following  set  of  transpositions 

{(l  2), (2  3),...,(nt -1  «<.)]•  For  G  the  size  of  each  of  these  neighborhoods  is 
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3.3.2  Inter-Orbit  Move  Neighborhoods 

Any  action  that  changes  the  orbit’s  characterization  matrix  causes  the  search 
to  leave  the  current  orbit  and  enter  a  new  one.  We  can  increment  or  decrement  the 
cardinality  of  one  or  more  of  the  non-zero  values  in  one  or  more  clusters.  What  we 
cannot  do  is  change  the  number  of  clusters  or  what  variables  are  assigned  to  each 
cluster.  Doing  so  would  change  the  partitioning  structure  as  opposed  to  which 
partition  we  are  currently  exploring. 
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3.3.3  Bounding  the  Solutions  in  an  Orbit 

We  can  use  the  objective  function  coefficients  and  the  orbit’s  characteristic 
matrix  to  calculate  an  upper  and  lower  bound  on  the  solutions  (not  necessarily 
feasible)  in  the  orbit.  Using  these  bounds,  it  may  be  possible  to  determine  if  an  orbit 
contains  any  feasible  solutions  without  actually  exploring  the  orbit.  We  may  also  be 
able  to  avoid  searching  an  orbit  if  its  solution  bounds  are  dominated  by  the  best 
solution  found  so  far. 

Let  Cjk  be  the  cost  associated  with  variable  j  in  cluster  k.  Let  nzk  be  the 
number  of  non-zero  elements  in  cluster  k  (i.e.  the  sum  of  column  k  in  O1’).  We  can 
calculate  the  upper  bound  objective  value  for  each  cluster  k  by  performing  the 
following: 

1.  Sort  (within  cluster  k )  the  Cjk  in  ascending  order 

2.  Set  j  =  nk-  nzk  +  1,  /  =  1,  and  zub*  =  0 

3.  While  Opik  *0 

3a.  zub*  =  £ub*  +  i*cjk 

3b.  7  =  7  +  1 

3c.  Cf ik  =  Cf  ik  -  1 

4.  i  =  i  +1 

5.  If  i  <  d,  go  to  3. 

The  upper  bound  for  the  orbit  is  ^=i  zmk  ■  We  can  calculate  the  lower  bound  of  an 
orbit  in  the  same  fashion  by  traversing  the  columns  of  Cf  from  bottom  to  top. 
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Example  3.2 

Assume  we  have  the  following  orbit  matrix 

|2  0  2l 
0  1  0 
0  0  1 
0  0  0 
op  =  0  1  0 
0  0  0 
0  0  0 
0  0  0 
[0  0  oj 

with  n  =  10,  ti\  =  3,  «2  =  3,  H3  =  4,  d  =10  and  r  =  3.  Let  our  cost  (sorted  by  cluster 
then  in  ascending  order)  be 

[l  2  4|  1  3  5|  2  2  4  6] 

The  upper  bound  for  cluster  1  is 

1.  zub1  =1*2  =  2 

2.  zub1  =  zub1  +  1  *4  =  6 

We  also  have  zm  =  31  and  zub3  =  24  so  zub  =  61 .  No  solution  in  the  orbit  feasible 
or  infeasible  will  have  a  solution  greater  than  61.  We  also  have  zlb  =  3  +  1 1  +  12  = 
26. 

We  may  be  able  to  get  tighter  bounds  by  representing  the  orbit  characteristics 
as  additional  constraints  and  solving  the  LP  relaxation  for  that  orbit.  Let  xjk  be  the  /th 
variable  in  cluster  k.  The  constraint  or  cut  associated  with  cluster  k  is 
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Adding  the  above  cut  for  each  cluster  and  solving  the  LP  relaxation  gives  us  an  upper 
or  lower  bound  depending  on  whether  we  are  maximizing  or  minimizing.  If  the  LP 
relaxation  is  infeasible  the  orbit  is  infeasible. 

If  we  are  maximizing  and  have  an  upper  bound  on  our  solution  from  some 
relaxation  method  and  the  lower  bound  of  the  orbit  is  greater  than  that  upper  bound, 
the  entire  orbit  is  infeasible  and  can  be  discarded.  If  we  have  a  current  best  solution 
greater  than  the  upper  bound  of  the  orbit  then  the  orbit  is  dominated  and  can  be 
discarded.  Similar  arguments  can  be  made  for  minimization. 

If  we  have  clustered  the  variables  based  on  reduced  costs  we  may  be  able  to 
put  a  bound  on  an  orbit  using  the  reduced  costs.  By  assigning  the  non-zero  values  in 
a  cluster  to  the  smallest  reduced  costs  of  the  variables  in  the  cluster  and  adding  across 
all  clusters,  we  can  calculate  the  minimum  weighted  distance  from  the  LP  optimal 
solution.  We  subtract  this  value  from  the  value  of  the  LP  optimal  solution  and  we 
have  an  upper  (max)  or  lower  (min)  bound  on  solutions  in  the  orbit. 

3.4  Orbit  Size 

Not  all  orbits  will  be  of  the  same  size.  The  size  of  the  orbit  is  based  on  the 
number  of  clusters  and  the  number  of  combinations  of  the  non-zero  elements  in  the 
clusters.  The  number  of  iterations  the  search  spends  in  an  orbit  should  be  a  function 
of  orbit  size. 

Again  let  d  be  the  order  of  the  vector  elements,  r  the  number  of  clusters  and  n* 
the  number  of  variables  assigned  to  cluster  k.  For  orbit  p  cluster  k,  the  number  of 
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ways  to  assign  the  first  non-zero  value  is  .  After  assigning  the  first  non-zero 

u) 


f  ~Opn  ) 

value  the  number  of  ways  to  assign  the  second  is  and  so  on.  So  the 

O 


number  of  combinations  for  orbit  p  cluster  k  is 


The  number  of  solutions  in  orbit  k  is  the  product  of  the  number  of  combinations  for 
each  cluster  in  the  orbit.  So  the  size  of  the  orbit  is 


PI  PI  ^  Hm=l0P> 


*=T  w  Op:„ 


Example  3.3 

For  the  following  orbit 


f2  0  2l 
0  1  0 
0  0  1 
0  0  0 
Op  =  0  1  0 
0  0  0 
0  0  0 
0  0  0 
[0  0  oj 


f  3  )  f3V2i  (4](2) 
the  size  is  •  •  =216. 

.2;  v  1 J  v  1  J  1  2  Jl  1, 
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3.5  Orbit  Landscapes 

As  stated  in  Section  2.2.3,  if  a  landscape  is  elementary  all  local  optima  are  at 
least  better  than  the  average  solution.  The  orbits  partitioning  the  solution  space  are 
distinct.  While  they  may  share  the  objective  function  and  neighborhood  definition, 
each  orbit  has  its  own  set  of  solutions.  A  neighborhood  that  creates  an  elementary 
landscape  in  each  orbit  implies  all  of  the  orbit’s  local  optima  are  better  than  the 
orbit’s  average  solution.  This  can  greatly  enhance  the  search  if  the  solution  space  is 
partitioned  into  “good”  and  “bad”  orbits. 

3.6  Escaping  Attractors 

Orbits  also  provide  a  deterministic  method  of  escaping  chaotic  attractors 
(Colletti  1999).  Battiti  and  Tecchiolli  (1994)  describe  chaotic  attractors  as  landscape 
characteristics  that  keep  the  search  confined  to  a  limited  region  of  the  solution  space. 
Reactive  tabu  search  is  designed  to  detect  and  escape  from  such  regions.  Such  escape 
procedures  can  cause  a  significant  change  in  the  current  solution.  Search  methods 
like  those  described  above  can  escape  the  attractor  simply  by  moving  to  a  different 
orbit. 


3.7  Parallel  Exploration 

Since  the  orbits  are  disjoint,  this  type  of  partitioning  lends  itself  well  to 
parallel  search.  A  centralized  control  algorithm  could  generate  orbit  characteristic 
matrices  and  pass  them  to  available  parallel  processors.  After  the  search  of  an  orbit  is 
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completed,  based  on  stopping  criteria,  the  processor  could  return  the  results  to  the 
control  algorithm  and  be  assigned  a  new  orbit. 

3.8  Conclusions 

Partitioning  the  solution  space  into  orbits  allows  the  algorithm  to  intensify  the 
search  in  areas  of  the  solution  space  believed  to  contain  good  solutions.  The  orbits 
keep  the  search  contained  in  these  areas  and  the  clusters  work  as  an  enhanced 
candidate  list,  reducing  the  total  number  of  moves  in  the  neighborhood  while  still 
retaining  the  “good”  moves.  In  the  next  chapter,  we  present  a  reactive  tabu  search 
implementation  of  these  concepts  for  the  unicost  set  covering  problem. 
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Chapter  4  -  A  Group  Theoretic  Tabu  Search  Algorithm  for 


Unicost  Set  Covering  Problems  (GTTS-USCP) 

In  this  chapter  we  develop  a  group  theoretic  tabu  search  (GTTS)  algorithm  for 
solving  the  unicost  set  covering  problem  (SCP),  the  GTTS-USCP.  We  solve  a  linear 
programming  (LP)  relaxation  of  the  problem  and  use  the  LP  optimum  to  construct  a 
quality  solution  profile.  As  described  in  the  previous  chapter,  we  use  group  theory  to 
partition  the  solution  space  into  orbits  based  on  this  profile.  We  tested  our  algorithm 
on  65  benchmark  problems  and  compared  the  results  against  the  previous  best  known 
and  solutions  obtained  by  CPLEX  9.0.  GTTS-USCP  discovered  46  new  best  known 
solutions.  GTTS-USCP  converged  significantly  faster  than  CPLEX  for  all  problem 
sets. 


4.1  Problem  Definition  and  Historical  Background 

The  set  covering  problem  (SCP)  is  a  well-known  combinatorial  optimization 
problem.  Given  a  0-1  incidence  matrix  A  with  m  rows  and  n  columns,  the  problem  is 
to  select  the  minimum  weight  subset  of  columns  while  ensuring  every  row  is  covered. 
Formally,  the  problem  is: 

Minimize  z  =  V"  wjc. 

j=  1  J  J 

Subject  to  X”=i  aijxi  - 1  i  =  1, ...,  m  (PI) 

Xj  e {0,1}  j  =  !,...,« 
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If  fly  =  1,  column  j  covers  row  i  and  wj  is  the  weight  or  cost  of  column  j,  ay. .  If  a;.  is 

selected  to  be  in  the  subset,  xj  is  set  to  1  and  wj  is  added  to  the  cost  of  the  solution. 
When  Wj  =  1  for  all  j,  the  problem  is  the  unicost  SCP  (Grossman  and  Wool  1997). 
The  cardinality  of  any  PI  solution,  *,  is  the  number  of  jc.=1,  which,  in  the  case  of  the 
unicost  SCP,  is  the  value  of  the  objective  function,  z.  The  linear  programming  (LP) 
relaxation  of  PI  where  the  x .  can  be  non-integer  is  denoted  by  PI  which  has  optimal 

solution,  x*p ,  and  optimal  objective  function,  z*p  . 

The  SCP  has  many  practical  applications  including  crew  scheduling  (Balas 
and  Carrera  1996,  Ceria  et  al.  1998,  Combs  2002,  Combs  and  Moore  2004), 
emergency  facility  location  (Daskin  and  Stem  1981,  Toregas  et  al.  1971)  and  political 
redistricting  (Garfinkel  1970).  The  SCP  is  known  to  be  NP-Hard  and  both  exact 
(Balas  and  Carrera  1996,  Balas  and  Ho  1980,  Beasley  1987,  Daskin  and  Stem  1981, 
Garfinkel  1970,  Toregas  et  al.  1971)  and  heuristic  (Beasley  1987,  Beasley  and  Chu 
1996,  Ceria  et  al.  1998,  Chvatal  1979,  Grossman  and  Wool  1997)  approaches  have 
been  proposed  for  it.  Several  of  these  have  made  extensive  use  of  Lagrangian 
relaxation  (Balas  and  Carrera  1996,  Balas  and  Ho  1980,  Beasley  1987,  Beasley  and 
Chu  1996,  Ceria  et  al.  1998)  and  column  dominance;  neither  of  which  are  as  effective 
for  unicost  SCPs. 

Grossmann  and  Wool  (1997)  explore  the  performance  of  nine  different 
heuristic  algorithms  on  the  unicost  SCP.  Most  algorithms  tested  were  LP  rounding- 
based  and  construction-based  approaches.  A  neural  network  algorithm  was  also 
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included  in  the  study.  A  randomized  version  of  the  greedy  construction  algorithm 
(Chvatal  1997,  Johnson  1974)  produced  the  best  results. 

4.2  Methodology 

4.2.1  Overview 

We  propose  to  solve  the  unicost  SCP  with  a  group  theoretic  tabu  search 
(GTTS-USCP)  algorithm.  Local  search  heuristics  provide  a  distinct  advantage  in 
solving  the  unicost  SCP  as  they  allow  the  search  to  be  based  on  the  cardinality  of  the 
solution.  As  described  in  chapter  3,  we  use  group  theory  to  partition  the  solution 
space  into  orbits  based  on  a  solution  profile  created  from  the  optimal  solution  to  the 
LP  relaxation.  The  LP  relaxation  of  an  IP  is  used  extensively  in  exact  methods,  but 
with  the  exception  of  the  LP  bound,  is  not  commonly  used  in  metaheuristics.  For  the 

SCP,  the  PI  solution  provides  valuable  information  for  a  direct  search  approach. 
The  choice  of  the  non-basic  variables,  xN ,  the  basic  variables,  xB ,  and  the  reduced 

costs  of  thex^,  cN,  can  provide  valuable  insight  into  the  characteristics  of  high- 
quality  IP  solutions. 

The  methodology  is  presented  in  detail  in  the  sections  that  follow.  However, 
for  the  purposes  of  clarity  and  ease  of  understanding,  we  provide  a  global  overview 

here.  First,  CPLEX  solves  PI .  The  associated  optimal  basic  variables,  the  xBi,  and 

the  cN  are  used  to  assign  the  Xj  to  clusters.  A  TS  starting  solution ,  x° ,  is  generated 
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by  setting  xB  =  0  and  then  selecting  a  subset  of  the  xBi  to  raise  to  value  1  so  that  all 
rows  are  covered  in  PI.  A  RTS  procedure  is  used  to  explore  the  orbits  with 
cardinality  between  that  of  x°  and  **p .  A  separate  RTS  procedure  is  also  used  to  find 
quality  orbits  to  explore.  The  algorithm  terminates  when  a  GTTS-USCP  solution 
value  equals  the  PI  lower  bound,  zLB  =  |"zPp ]  or  when  the  allotted  time  has  expired. 

4.2.2  Clustering  Variables  Based  on  the  LP  Relaxation 

The  CPLEX  solution  to  PI  details  z*L P,  xB,  ~cN ,  and  which  xN  are  at  their 
upper  (1)  or  lower  (0)  bounds.  All  xBi  join  the  same  cluster.  The  upper  bound 
variables,  xNU] ,  and  lower  bound  variables,  xNLj,  cluster  separately  based  on  their 

reduced  costs,  cNUJ  and  cNLj.  For  unicost  PI ,  all  cNLj  <1.  If  cNLj  =1,  ay.  does  not 
cover  any  of  rows  associated  with  the  binding  constraints  at  optimality.  A  cNUj  =  -1 
implies  that  if  aj  is  removed  from  the  solution,  we  must  select  at  least  two  new 
columns  to  render  the  LP  solution  feasible.  While  it  is  possible  for  cNUj  to  be  less 
than  -1,  it  is  quite  unusual.  To  avoid  creating  too  many  clusters  and  too  many  orbits, 
we  place  the  xNJ  in  the  same  cluster  if  their  cNj  are  equal  to  the  nearest  0.1  digit.  As 

illustrated  in  Figure  4.1,  the  clusters  are  created  in  the  following  order:  upper  bound 
variables,  basic  variables,  lower  bound  variables. 
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Upper  bound 

Basic 

Lower  bound 

...  ...  ... 

-0.2 

-0.1 

-0.0 

0.0 

0.0 

0.1 

0.2 

1.0 


Figure  4.1  -  Variable 


We  initially  set  all  xNUj  =  1 ,  all  xNLj  -  0  and  all  xBj  =  0  making  the  solution 

integer  feasible.  Integer  feasibility  is  then  maintained  throughout  the  algorithm. 
Since  the  problem  is  binary,  the  elements  of  the  solution  are  order  2.  The  orbit 
characterization  matrix  is  then  a  1  by  r  vector,  where  r  is  the  number  of  clusters.  For 
example,  a  solution  to  a  22  variable  problem  could  be 

[l  1 1  l|l  101|10100|001  [oooooo] 

where  the  solution  cardinality  is  z  =  10  and  the  orbit  characterization  vector  is 
Op  =(4  3  2  1  0). 


4.2.3  Partitioning  the  Solution  Space  into  Orbits 

The  solution  space  is  partitioned  into  orbits  first  by  total  solution  cardinality, 
z,  and  second  by  the  cardinality  of  the  individual  clusters  as  described  in  Op.  Initially 
all  orbits  are  open.  A  list  of  orbits  visited  is  maintained  and  an  orbit  is  closed  to 
further  search  after  it  has  been  searched  for  MAX_ORBIT_ITER  iterations  or  has 
been  visited  MAX_ORBIT_VISITS  times.  An  orbit  hash  function  (Woodruff  and 
Zemel  1993),  Q.(k) ,  is  used  to  facilitate  access  to  orbit  information.  Orbit  hash 
values  are  calculated  by  generating  a  random  number,  x  ■ ,  for  each  cluster  j.  The 

orbit’s  hash  value  is  Q (Op )  =  ^.=1  TjOPj  .  Orbits  can  be  further  distinguished  by  the 
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number  of  free  ones  and  free  zeros  in  the  orbit.  Free  ones  (zeroes)  are  the  number  of 
ones  (zeroes)  in  clusters  that  are  not  all  ones  (zeroes).  For  example,  the  solution 
below  has  6  free  ones  and  6  free  zeroes. 

[l  1 1 1  |l  1 0 1  |l 0 1 00 100 1  |oooooo] 

Any  orbit  with  Op  j  less  than  zLB  contains  only  infeasible  solutions  and 

may  be  immediately  closed.  Once  a  feasible  solution  is  found  with  cardinality  zUB , 
all  orbits  with  0Ps  >  zUB  are  dominated  and  may  be  closed.  It  is  probable  that 

orbits  with  large  Op j  in  the  upper  bound  clusters  and  small  Op.  in  the  lower  bound 

clusters  will  contain  good  solutions.  This  partitioning  scheme  permits  the 
concentration  of  search  effort  within  these  good  orbits. 

4.2.4  Inter-  and  Intra-Orbit  Neighborhoods 

As  shown  in  Section  3.3.1,  our  group  acting  on  the  set  of  binary  vectors  size  n 

is  G  =  Snlx...xSnr  where  , nr  =  n .  For  our  intra-orbit  neighborhood  we  need  a 

set  H  c  G  such  that  G  =  .  We  let  H  equal  the  set  of  all  transpositions  in  G.  This 

is  equivalent  to  all  moves  swapping  two  elements  in  the  same  cluster. 

Any  move  that  changes  Op  is  an  inter-orbit  move.  We  use  3  different  inter¬ 
orbit  neighborhoods  at  different  points  in  the  GTTS-USCP  algorithm.  Op .  is 
increased  {select  neighborhood)  or  decreased  {unselect  neighborhood)  by  “toggling” 
the  value  of  a  single  cluster  j  variable.  Both  Op]  and  Opk  may  be  changed  by 
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swapping  values  of  variables  in  clusters  j  and  k.  When  searching  inter-orbit 
neighborhoods  we  consider  only  open  orbits. 

For  all  neighborhoods,  the  best  move  is  defined  in  terms  of  deficit,  the  number 
of  unsatisfied  rows,  and  surplus,  the  total  amount  that  rows  are  over-satisfied.  A  row 
is  over-satisfied  when  it  is  covered  more  than  once.  A  row  covered  by  3  columns  has 
a  surplus  of  2.  Once  a  feasible  solution  is  found  we  decrease  z  by  unselecting  a 
column,  therefore  GTTS-USCP  is  usually  searching  an  orbit  for  feasible  solutions. 
The  best  move  yields  the  smallest  deficit.  Deficit  ties  are  broken  by  largest  surplus 
and  surplus  ties  are  broken  by  order  of  evaluation  (lexicographically). 

A  complete  swap  neighborhood  (which  ignores  cluster  membership)  would  be 
0(  n2 ).  For  intra-orbit  swaps,  this  effort  is  reduced  by  considering  swap  pairs  only  if 
they  are  in  the  same  cluster.  Further,  if  Op ;  =  nj  or  Of  =  0  no  swaps  are  considered 

in  cluster  j.  However,  even  with  this  reduction  a  complete  within-cluster  or  outer- 
cluster  swap  neighborhood  is  too  costly  for  large  problems.  Both  the  intra-orbit  and 
inter-orbit  swap  neighborhoods  are  based  on  a  conditional  select  move.  First,  moves 
that  unselect  a  non-tabu  selected  column  in  the  current  solution  are  evaluated.  The 
best  of  such  moves  is  chosen.  Given  that  chosen  column  will  be  unselected,  we  next 
find  the  best  non-tabu  unselected  column  to  select. 

A  complete  select  neighborhood  is  O(n).  This  neighborhood’s  size  may  be 
reduced  by  incorporating  a  candidate  list  heuristic  that  also  helps  diversify  the  search. 
A  column  selection  is  considered  only  when  the  current  solution  is  infeasible.  We 
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first  find  the  row  that  has  been  unsatisfied  for  the  longest  number  of  iterations,  then 
we  only  select  from  the  columns  that  will  satisfy  that  row. 

Since  ties  in  surplus  are  broken  lexicographically,  the  order  in  which  the 
clusters  are  evaluated  affects  the  performance  of  the  algorithm.  For  select 
neighborhoods,  we  first  examine  the  upper  bound  clusters  then  the  basic  cluster  and 
finally  the  lower  bound  clusters.  This  favors  increasing  the  cardinality  of  the  upper 
bound  clusters  over  the  rest.  Unselect  neighborhoods  are  examined  in  the  opposite 
order  of  select  neighborhoods  favoring  decreasing  the  cardinality  of  the  lower  bound 
clusters  over  the  rest. 

4.2.5  Tabu  Lists  and  Tabu  Tenure 

Two  types  of  tabu  structures  are  used  in  our  algorithm.  The  first  is  a  broad- 
gauge  structure  which  tracks  the  last  iteration  a  column  was  selected  or  unselected.  A 
selected  or  unselected  column’s  status  cannot  be  changed  again  for  tabu  tenure 
iterations.  The  tabu  tenure  may  increase  or  decrease  when  the  algorithm  detects 
cycling. 

The  second  type  of  tabu  structure  is  a  fine-gauge  structure  which  tracks  each 
individual  solution,  noting  when  it  was  last  visited,  and  how  many  times  it  has  been 
visited.  For  efficiency,  another  hash  function  (Woodruff  and  Zemel  1993),  cp(x),  is 
used.  A  random  number,  ps ,  is  generated  for  each  Xj .  The  solution’s  hash  value  is 
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<P(x)~’m]-iPjxj  ■  Solutions  can  be  further  distinguished  by  their  deficit  and 
surplus. 

Each  orbit  maintains  its  own  tabu  structures.  In  addition,  a  separate  tabu 
structure  is  used  during  the  RTS  procedure  to  find  quality  orbits.  The  default  tabu 
tenure  for  select  moves  is  SELECT_TENURE  and  the  default  tabu  tenure  for  unselect 
moves  is  UN SELECT_TENURE . 

The  solution  tabu  structure  is  used  to  detect  cycling  and  to  control  the  tabu 
tenures.  If  a  solution  is  repeated,  the  tabu  tenures  are  increased  by  a  multiplicative 
factor  (*1.618).  If  MIN_NEW_SOLS  consecutive  new  solutions  are  visited,  the  tabu 
tenures  are  returned  to  their  default  values.  If  MAX_REPEATED_SOLS  solutions 
are  repeated  MAX_REPEATS  times,  the  search  is  presumed  to  be  in  a  chaotic 
attractor  basin  (Battiti  and  Tecchiolli  1994)  and  an  escape  is  achieved  by  departing 
the  current  orbit  or  by  increasing  the  tabu  tenure  if  we  are  still  searching  for  an  orbit. 

4.2.6  Finding  a  Starting  Solution 

To  obtain  a  starting  solution,  all  xNj  are  fixed  at  their  values  in  jt*p  and  the 
xBj  are  set  to  zero.  As  illustrated  in  the  pseudo  code  of  Figure  4.2,  basic  columns  are 
then  selected  until  zLB  is  reached.  If  that  solution  is  not  feasible,  intra-orbit  swap 
moves  are  performed  within  the  basic  cluster,  until  no  improving  move  is  available 
(steepest  descent).  If  the  solution  is  still  not  feasible,  a  select  move  is  performed 
within  the  basic  cluster  and  the  process  repeats  until  feasibility  is  achieved.  Since 
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selecting  all  columns  would  certainly  be  a  feasible,  albeit  poor,  solution  we  will 

converge  to  a  feasible  solution  in  a  maximum  of  n  -  zLB  steps.  After  feasibility  is 

achieved,  any  identified  redundant  columns,  columns  whose  unselection  will  not 

destroy  feasibility,  are  removed. 

Generate  Feasible  Starting  Solution 

Select  columns  from  the  basic  cluster  until  the  zLB  is  reached 
If  feasible 

Terminate  with  optimal  solution 

Else  { 

While  not  feasible  { 

Select  a  column  from  the  basic  cluster 
While  improving  move  found 

Execute  intra-orbit  swap  move 

} 

Attempt  to  unselect  redundant  columns 

} 

Save  the  solution  and  zUB 

Figure  4.2  -  Starting  Solution  Algorithm 


4.2.7  Primary  Search  Strategy 

A  pseudo-code  of  the  primary  search  strategy  is  presented  in  Figure  4.3. 
After  the  initial  solution  is  obtained,  the  GTTS-USCP  still  focuses  on  the  basic 
variables.  A  non-basic  variable’s  status  is  modified  only  if  doing  so  achieves 

feasibility.  Since  the  duality  gap,  zUB-  zLB,  is  known,  [~(zUB  _^lb)/2]  columns  are 

unselected  from  the  basic  cluster.  Next,  an  RTS  procedure,  based  on  the  intra-orbit 
swap  neighborhood,  is  applied  to  the  resulting  orbit  until  (1)  a  feasible  solution  is 
found,  (2)  MAX_ORBIT_ITER  iterations  have  been  performed,  (3) 


MAX_NI_ORBIT_ITER  iterations  have  been  performed  without  improving  the  best 
orbit  solution,  or  (4)  the  time  limit  is  reached.  If  a  feasible  solution  is  obtained,  the 

duality  gap  is  recalculated,  [~(zUB  -  zLB  )/2]  columns  are  unselected  from  the  basic 
cluster  and  the  process  repeats. 

If  a  feasible  solution  is  not  found,  the  search  departs  the  current  orbit  and  the 
orbit’s  best  solution  is  instantiated  as  the  current  incumbent  solution.  Next,  the 
current  inter-orbit  swap  neighborhood  is  evaluated  in  pursuit  of  a  move  leading  to  a 
feasible  solution.  If  that  search  yields  feasibility,  the  duality  gap  is  recalculated  and 
the  process  repeats;  if  not,  the  current  select  neighborhood  is  evaluated.  If  feasibility 
is  achieved,  the  duality  gap  is  recalculated  and  the  process  repeats.  If  both 
neighborhoods  fail  to  find  a  feasible  solution,  a  basic  cluster  column  is  selected  and 
the  new  orbit  is  searched.  If  zUB  is  reached  without  achieving  feasibility,  the  search 
is  expanded. 
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Primary  Search 

Search  Orbit  with  a  RTS  procedure 

Perform  intra-orbit  swap  moves  until  a  feasible  solution  is  found  or 
termination  criteria  is  met 
If  feasible  { 

Attempt  to  unselect  redundant  columns 
Save  the  solution  and  zUB 

If  ~UB  —  Ab 

Terminate  with  optimal  solution 

Else  { 

Remove  ( zUB  -  zLB )  /  2  columns  from  the  basic  cluster 
Goto  Primary  Search 

} 

} 

Else  { 

If  feasible  inter-orbit  swap  move  found  { 

Execute  swap 
Goto  Primary  Search 

} 

Else  If  z  + 1  <  zUB  { 

If  feasible  select  column  move  found  { 

Select  the  column 
Goto  Primary  Search 

} 

Else  { 

Select  a  column  from  the  basic  cluster 
Goto  Primary  Search 

} 

} 

Else 

Goto  Expanded  Search 


Figure  4.3  -  Basic  Search  Algorithm 
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Expanded  Search 

Search  for  a  feasible  solution  with  a  RTS  procedure 

Perform  intra-orbit  and  inter-orbit  swap  moves  until  a  feasible  solution 
is  found  or  termination  criteria  is  met 
If  feasible  { 

Attempt  to  unselect  redundant  columns 
Save  the  solution  and  zm 

If  £ub  =  Ab 

Terminate  with  optimal  solution 

Else  { 

Unselect  a  column 

While  feasible  solution  found  { 

Search  Orbit  with  Reactive  Tabu  Search 

Perform  intra-orbit  swap  moves  until  a  feasible 
solution  is  found  or  termination  criteria  is  met 
If  feasible  { 

Attempt  to  unselect  redundant  columns 
Save  the  solution  and  zUB 

If  ^UB  —  ^LB 

Terminate  with  optimal  solution 
Unselect  a  column 

} 

} 

Goto  Expanded  Search 

} 

} _ 


Figure  4.4  -  Expanded  Search  Algorithm 


4.2.8  Expanding  the  Search 

A  pseudo-code  of  the  expanded  search  strategy  is  presented  in  Figure  4.4. 
After  exploring  the  orbits  near  the  optimal  LP  solution,  GTTS-USCP  expands  the 
search  to  other  areas  of  the  solution  space.  A  RTS  procedure,  based  on  both  inter¬ 
orbit  and  intra-orbit  swap  neighborhoods,  is  used  to  find  a  good  region  for 
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exploration.  This  RTS  procedure  continues  until  either  a  feasible  solution  is  found  or 
the  time  limit  is  reached. 

If  a  feasible  solution  is  found,  we  unselect  a  single  column  (from  any  cluster) 
and  explore  the  resulting  orbit.  If  a  feasible  solution  is  found  while  exploring  the 
orbit,  we  unselect  a  column  again  and  the  process  repeats.  If  a  feasible  solution  is  not 
found  while  exploring  the  orbit,  we  begin  another  RTS  procedure  at  the  new  z-  The 
process  ends  when  the  time  limit  is  reached  or  a  solution  equal  to  zLB  is  found. 

4.3  Computational  Results 

4.3.1  Test  Cases 

The  benchmark  problems  solved  were  obtained  from  Beasley’s  OR-Library 
(Beasley  1990).  Problem  sets  4-6  originally  appeared  in  (Balas  and  Ho  1980), 
problem  sets  A-D  appeared  in  (Beasley  1987)  and  problem  sets  NRE-NRH  appeared 
in  (Beasley  1990).  All  problems  were  randomly  generated  based  on  the  strategy  of 
(Balas  and  Ho  1980).  All  of  these  problems  were  generated  as  weighted  SCPs. 
Grossman  and  Wool  (1997)  solved  all  but  NRG  and  NRH  as  uni  cost  SCPs.  One 
problem  set,  E  from  (Beasley  1987),  is  unicost  but  is  not  solved  here  due  to  its  trivial 
size  (all  algorithms  reach  an  optimal  solution  in  0  seconds). 
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4.3.2  Test  Procedures 


All  tests  were  performed  on  Dell  Precision  530  Workstations  running  SuSE 
Linux  with  two  1.8GHz  Pentium  Xeon  processors  and  1GB  of  RAM.  The  machines 
are  multi-user  platforms.  An  attempt  was  made  to  find  machines  that  were  not  too 
busy,  but  as  each  problem  ran  for  at  least  an  hour,  CPU  usage  surely  fluctuated  during 
processing.  Each  problem  was  solved  using  CPLEX  9.0  and  GTTS-USCP.  The  time 
limits  used  were  7200  seconds  for  CPLEX  and  3600  seconds  for  GTTS-USCP.  The 
GTTS-USCP  algorithm  was  coded  in  C. 

The  previously  published  best  known  solutions  for  these  problems  were 
published  in  Grossman  and  Wool  (1997).  They  performed  their  tests  on  an  IBM 
RS6000  model  370  workstation  with  128MB  of  RAM.  They  also  coded  their 
algorithms  in  C. 

4.3.3  CPLEX  9.0 

The  algorithms  tested  by  Grossman  and  Wool  (1997)  are  unsophisticated  by 
today’s  standards.  To  provide  a  more  modern  benchmark,  we  compare  our  results  to 
CPLEX  version  9.0.  CPLEX  uses  a  very  sophisticated  branch  and  cut  algorithm  to 
solve  mixed  integer  programs  (MILPs).  The  algorithm  is  further  aided  by  two 
heuristics.  The  first  attempts  to  create  a  feasible  solution  from  the  fractional  solution 
at  the  node.  The  second  attempts  to  improve  the  incumbent  integer  solution  through  a 
neighborhood  search  (ILOG  2003). 
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CPLEX  was  also  used  to  solve  PI  for  the  GTTS-USCP  algorithm.  The  dual 
simplex  LP  solver  was  used  for  the  smaller  problem  sets,  4-6  and  A-D.  The  sifting 
LP  solver  was  used  for  the  larger  problem  sets,  NRE-NRH.  The  default  settings  were 
used  for  all  other  parameters. 

4.3.4  Results 

Tables  4.1  and  4.2  contain  the  results  of  our  tests  as  well  as  the  problem 
details  and  previous  best  known  solutions.  The  best  solution  found  for  each  problem 
is  highlighted  in  bold.  GTTS-USCP  found  the  best  solution  on  59  of  the  65 
problems.  It  outperformed  CPLEX  on  47  of  65.  It  significantly  outperformed 
CPLEX  on  the  larger  problem  sets  NRG  and  NRH.  Curiously,  neither  CPLEX  nor 
GTTS-USCP  were  able  to  do  as  well  as  R-Gr  (Grossman  and  Wool  1997)  on  the 
problem  sets  with  higher  density,  NRE  and  NRF.  None  of  the  solutions  were  proven 
to  be  optimal. 

It  should  be  noted  here  that  CPLEX  was  executed  using  the  default  settings 
while  GTTS-USCP,  as  well  as  the  algorithms  from  Grossman  and  Wool  (1997),  were 
specifically  designed  to  solve  the  unicost  SCP.  A  researcher  well-versed  in  CPLEX’ s 
parameters  and  settings  would  likely  be  able  to  improve  its  performance  through 
experimentation.  However,  the  disparity  in  performance  between  CPLEX  and  GTTS- 
USCP  is  quite  dramatic  and  we  do  not  believe  any  such  improvements  would  be 
enough  to  close  this  gap. 
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problem 

number 

of 

rows 

number 

of 

columns 

density 

previous 
best  known 

CPLEX 9 
best 

(seconds) 

GTTS- 

USCP 

best 

(seconds) 

4.1 

200 

1000 

2% 

41 

38(11) 

38  (938) 

4.2 

200 

1000 

2% 

38 

37  (42) 

37(5) 

4.3 

200 

1000 

2% 

40* 

38  (28) 

38(1) 

4.4 

200 

1000 

2% 

41 

39(851) 

38  (272) 

4.5 

200 

1000 

2% 

40 

39  (64) 

38  (23) 

4.6 

200 

1000 

2% 

40 

38  (50) 

32(3) 

4.7 

200 

1000 

2% 

41 

39(211) 

38(413) 

4.8 

200 

1000 

2% 

40 

38(131) 

38(6) 

4.9 

200 

1000 

2% 

40 

38  (835) 

38(35) 

4.10 

200 

1000 

2% 

41 

38(1772) 

38(161) 

problem  set  4  average 

40.2 

38.2  (399.5) 

37.8  (185.7) 

5.1 

200 

2000 

2% 

35 

35  (824) 

35(5) 

5.2 

200 

2000 

2% 

35 

35(1.37) 

35(6) 

5.3 

200 

2000 

2% 

36 

35  (373) 

34(39) 

5.4 

200 

2000 

2% 

36 

35  (122) 

34(1182) 

5.5  . 

200 

2000 

2% 

36 

35  (2257) 

34(12) 

5.6 

200 

2000 

2% 

36 

36  (29) 

34  (989) 

5.7 

200 

2000 

2% 

35* 

35  (33) 

34(75) 

5.8 

200 

2000 

2% 

37 

35  (85) 

34  (74) 

5.9 

200 

2000 

2% 

36 

36(113) 

35  (6) 

5.10 

200 

2000 

2% 

36 

35  (4739) 

34(1873) 

problem  set  5  average 

35.8 

35.2  (871.2) 

34.3  (426.1) 

6.1 

200 

1000 

5% 

21 

22  (3) 

21(5) 

6.2 

200 

1000 

5% 

21 

21(2) 

21(6) 

6.3 

200 

1000 

5% 

21* 

22  (3) 

21  (10) 

6.4 

200 

1000 

5% 

22 

22  (40) 

21(4) 

6.5 

200 

1000 

5% 

22 

22  (2) 

21  (25) 

problem  set  6  average 

21.4 

21.8(10) 

21  (10) 

Table  4.1  -  Results  for  problem  sets  4-6 
*  -  Solution  found  by  Neural  Network  algorithm  (Grossman  and  Wool  1997) 
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problem 

number 

of 

rows 

number 

of 

columns 

density 

previous 
best  known 

CPLEX 9 
best 

(seconds) 

GTTS-USCP 

best 

(seconds) 

A.l 

300 

3000 

2% 

40 

42  (19) 

39  (337) 

A.2 

300 

3000 

2% 

41 

41  (21) 

39(79) 

A.3 

300 

3000 

2% 

40 

40  (398) 

39  (179) 

A.4 

300 

3000 

2% 

40 

42  (19) 

38(1715) 

A.5 

-  300 

3000 

2% 

40 

39  (4598) 

38(771) 

problem  set  A  average 

40.2 

40.8(1011) 

38.6  (616.2) 

B.l 

300 

3000 

5% 

23 

23  (826) 

22(719) 

B.2 

300 

3000 

5% 

22 

23  (134) 

22(17) 

B.3 

300 

3000 

5% 

22 

23  (12) 

22  (698) 

B.4 

300 

3000 

5% 

23 

23  (77) 

22(1910) 

B.5 

300 

3000 

5% 

23 

23  (2422) 

22(46) 

problem  set  B  average 

22.6 

23  (694.2) 

22  (678) 

C.l 

,  400 

4000 

2% 

45 

48  (3251) 

43  (1524) 

C.2 

400 

4000 

2% 

45 

46  (4488) 

44(197) 

C.3 

400 

4000 

2% 

45 

49  (41) 

43  (1029) 

C.4  . 

'  .  400 

4000 

2% 

46 

47  (3994) 

43  (1325) 

C.5 

.  400 

4000 

2% 

45 

48  (6006) 

44(149) 

problem  set  C  average 

45.2 

47.6  (3556) 

43.4  (844.8) 

D.l 

400 

4000 

5% 

26 

27  (2823) 

25  (395) 

D.2 

400 

4000 

5% 

25A 

26  (32) 

25(1890) 

D.3 

400 

4000 

5% 

25 

26  (445) 

25(91) 

D.4 

400 

4000 

5% 

26 

26  (178) 

25  (226) 

D.5 

400 

4000 

5% 

26 

26  (3489) 

25  (200) 

problem  set  D  average 

25.8 

26.2(1393.4) 

25  (560.4) 

NRE.l 

500 

5000 

10% 

17 

18(77) 

18  (38) 

NRE.2 

500 

5000 

10% 

17 

18(77) 

18  (27) 

NRE.3 

,.  500 

5000 

10% 

17 

18  (455) 

18  (32) 

NRE.4 

500 

5000 

10% 

17 

18(71) 

12(54) 

NRE.5 

500 

5000 

10% 

17 

17  (586) 

18(176) 

problem  set  NRE  average 

17 

17.8  (255.2) 

17.8  (65.4) 

NRF.l 

■  500 

5000 

20% 

10 

1 1  (90) 

11  (29) 

NRF.2 

500 

5000 

20% 

11 

11  (94) 

11  (39) 

NRF.3 

500 

5000 

20% 

11 

11(132) 

11  (30) 

NRF.4 

,  500 

5000 

20% 

11 

11  (1083) 

11  (22) 

NRF.5 

.  500 

5000 

20% 

11 

10  (482) 

11  (23) 

problem  set  NRF  average 

10.8 

10.8  (376.2) 

1 1  (28.6) 

NRG.l 

1000 

10000 

2% 

- 

74  (582) 

63(1089) 

NRG.2 

1000 

10000 

2% 

- 

76(175) 

61(3401) 

NRG. 3 

1000 

10000 

2% 

- 

75  (6426) 

62  (901) 

NRG.4 

1000 

10000 

2% 

- 

74  (3723) 

63(1045) 

NRG.5 

1000 

10000 

2% 

- 

73  (577) 

63  (406) 

problem  set  NRG  average 

- 

74.4  (2296.6) 

62.4  (1368.4) 

NRH.l 

1000 

10000 

5% 

- 

40  (756) 

35  (2008) 

NRH.2 

1000 

10000 

5% 

- 

39  (5716) 

36  (297) 

NRH.3 

1000 

10000 

5% 

- 

39  (1527) 

36  (968) 

NRH.4 

1000 

10000 

5% 

- 

40  (748) 

35  (940) 

NRH.5 

1000 

10000 

5% 

- 

37  (5734) 

36  (454) 

problem  set  NRH  average 

- 

39  (2896.2) 

35.6  (933.4) 

Table  4.2  -  Results  for  problem  sets  A-D  and  NRE-NRH 
A  -  Solution  found  by  Alternating  Greedy  algorithm  (Grossman  and  Wool  1997) 
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Since  GTTS-USCP  and  CPLEX  were  executed  on  the  same  platform,  we  can 
make  an  accurate  comparison  of  the  convergence  properties.  For  the  purposes  of  this 
dissertation  we  define  convergence  as  the  time  or  effort  required  to  reach  a  similar 


level  of  solution  quality.  Table  4.3  shows  the  average  CPU  seconds  for  each  problem 
set  for  GTTS-USCP  and  CPLEX.  Since  the  algorithms  achieved  different  quality 


solutions,  the  time  to  the  best  common  solution  is  used  for  the  comparison  whenever 
possible.  When  no  common  quality  solution  exists,  the  next  lower  solution  for 
GTTS-USCP  is  used.  GTTS-USCP  converged  significantly  faster  than  CPLEX  for 
all  problem  sets. 

Table  4.4  compares  the  convergence  properties  of  GTTS-USCP  to  R-Gr 
(Grossman  and  Wool  1997).  R-Gr  was  executed  on  a  slower  computer  and  for  5  of 
the  problems  GTTS-USCP  did  not  reach  the  same  quality  solution.  As  expected  R-Gr 
is  faster  than  GTTS-USCP;  however,  GTTS-USCP  significantly  outperforms  R-Gr  in 
terms  of  overall  solution  quality.  If  R-Gr  was  executed  for  a  longer  period  of  time  or 


for  more  iterations,  it  is  unlikely  to  improve  upon  the  solutions  it  has  already  found. 


problem 

set 

CPLEX  9 
avg  (sec) 

GTTS-USCP 
avg  (sec) 

4 

38.2(116.6) 

5 

EHIB 

35.2  (4.4) 

6 

21.8(10) 

21.8  (2.8) 

A 

40.8(1011) 

40.8(21.2) 

B 

23  (694.2) 

23  (26.8) 

C 

47.6  (3556) 

47.4  (7.2) 

D 

26.2(1393.4) 

26  (46.2) 

NRE 

18  (150.6) 

18(57.8) 

NRF 

11  (296.2) 

11  (27) 

NRG 

74.4  (2296.6) 

68.6(114.4) 

NRH 

39  (2896.2) 

38.2(175) 

Table  4.3  -  CPLEX  vs  GTTS-USCP  solve  time 
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problem 

set 

R-Gr 

Avg  (sec) 

GTTS-USCP 
avg  (sec) 

4 

40.3  (1.6) 

39.6  (0.8) 

5 

35.9  (3) 

35.8  (2.6) 

6 

21.6  (2) 

21.6  (3.4) 

A 

40.2  (6) 

40  (7.8) 

B 

22.6  (8) 

22.6  (167.4) 

C 

45.2  (10) 

45.2  (20.8) 

D 

25.8  (14.2) 

25.8  (54.8) 

NRE 

17  (38) 

17.8  (65.4) 

NRF 

10.8  (71.2) 

11  (27) 

Table  4.4  -  R-Gr  vs  GTTS-USCP  solve  time 


4.4  Conclusion 

The  use  of  variable  clustering  and  group  theory  allowed  our  algorithm  to 
intensify  the  search  in  the  areas  of  the  solution  space  believed  to  contain  good 
solutions.  The  orbits  kept  the  search  contained  in  these  areas  and  the  clusters  worked 
as  an  enhanced  candidate  list,  reducing  the  total  number  of  moves  in  the 
neighborhood  while  still  retaining  the  “good”  moves.  These  techniques  proved  very 
effective  for  the  unicost  SCP  discovering  46  new  best  known  solutions  to  the 
benchmark  problems.  However,  these  techniques  are  very  problem  dependent.  In  the 
next  chapter  we  examine  group  theoretic  local  search  techniques  for  the  general  IP. 
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Chapter  5  -  Local  Search  Neighborhoods  for  General  IP 

The  feasible  region  for  a  linear  program  (LP)  is  a  convex  polyhedron  ( P ) 
which  may  be  bounded  or  unbounded.  If  the  LP  decision  variables  are  constrained  to 
be  integer,  the  result  is  an  integer  linear  programming  problem  (ILP).  The  convex 
hull  of  P,  conv(P),  is  the  smallest  polyhedron  containing  the  integer  points  of  P.  If 
we  can  define  conv(P),  we  can  relax  the  ILP  integrality  requirements  and  solve  the 
corresponding  LP.  Our  solution  will  be  integer  and  we  will  have  the  optimal  solution 
to  the  original  ILP.  Unfortunately  defining  conv(P)  is  impractical  for  most  problems. 

Often  we  are  able  to  find  a  feasible  integer  solution  within  P.  We  can  perform 
a  local  search  by  starting  from  this  solution  and  moving  to  an  adjacent  feasible  point 
by  increasing  or  decreasing  a  variable  value  and  checking  the  constraint  set  to  ensure 
we  remain  in  P. 

The  group  minimization  problem  (GMP)  developed  by  Gomory  (1965,  1967, 
1969)  is  the  mathematical  formulation  for  deriving  the  optimal  ILP  solution  from  the 
optimal  solution  to  the  corresponding  LP  (Johnson  1980).  Under  certain  conditions 
(detailed  below),  the  optimal  solution  to  the  GMP  is  the  optimal  solution  to  the 
original  ILP.  When  we  solve  the  GMP,  we  solve  the  ILP  and  vice-versa. 
Furthermore,  the  form  of  the  GMP  is  the  same  regardless  of  the  context  of  the 
original  ILP,  so  an  algorithm  to  solve  the  GMP  can  also  solve  any  ILP. 

The  feasible  region  for  the  GMP  is  also  a  convex  polyhedron,  called  the 
corner  polyhedron  (CP).  The  objective  function  costs  in  the  GMP  are  the  reduced 
costs  of  the  non-basic  variables  from  the  optimal  LP  solution.  These  costs  are  the 
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penalties  for  moving  away  from  the  LP  optimal  point,  x[p  along  each  non-basic 
variable’s  axis.  So  xpp  acts  as  an  anchor  point  for  our  search.  By  examining  the 
local  search  in  terms  of  the  GMP,  we  gain  insights  that  will  help  us  solve  the  ILP. 

5.1  The  Group  Minimization  Problem  (GMP) 

5.1.1  Derivation  of  the  GMP 

Given  an  ILP,  we  derive  the  associated  GMP  by  first  solving  the  LP 
relaxation.  Consider  the  ILP  (5.1). 

Max  cx 

s.t.  Ax  <  b  (5.1) 

x  >  0  and  integer 

where  c  is  a  size  n  row  vector,  jc  is  a  size  n  column  vector,  A  is  a  mxn  integer 
matrix  and  b  is  a  size  m  integer  column  vector.  Adding  positive  slack  variables  to 
change  the  constraints  to  equalities  yields  (5.2)  where  c  =  [c|  0],  x'  =  [x'|  s'],  and 

a  =  [A|/]. 

Max  cx 

s.t.  Ax  =  b  (5.2) 

x  >  0  and  integer 

Let  B  be  the  optimal  basis  from  the  LP  relaxation  and  N  be  the  matrix 
containing  the  non-basic  columns,  separating  A  into  B  and  N  yields  (5.3)  where 
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and  Xb  are  the  cost  coefficients  and  variables  associated  with  the  basic  columns  and 
and  x n  are  those  associated  with  the  non-basic  columns. 


Max  cBxB  +cNxN 
s.t.  Bxb  +  Nxn  =  b 

xB  >  0  and  integer 
xN  >  0  and  integer 


(5.3) 


Solving  for  xb  in  terms  of  xN  yields 


=  B  *  (b-NxN)  =  B  Ib-B1NxN 


(5.4) 


Substituting  for  Xb  using  Equation  5.4  yields 

Max  CgB-'b-^gB-'N -cN)xN 

s.t.  xB  -B'lb-B'J Nxn  ^ 

xB  >  0  and  integer 
xN  >  0  and  integer 

The  reduced  costs  from  the  optimal  LP  solution  are  defined  as 
cN  =  [cgB'1  N -cNy  Since  we  are  maximizing  and  the  variables  do  not  have  upper 

bounds,  cN  >  0  at  optimality  and  the  objective  function  value  is  z*LP  =  clsB'‘b . 
Therefore,  the  objective  function  in  (5.5)  can  be  rewritten  as  z  =  z*LP  ~cNxN  and  the 
problem  can  be  restated  as 

Min  cNxN 

s.t.  xR  =B'Ib-B'1  Nxn 

B  N  (5.6) 

xB  >  0  and  integer 

xN  >  0  and  integer 
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The  basic  variables,  xB,  will  be  integer  if  they  satisfy  the  congruence 
relationship  B^b -B'1  Nxn  =0  (mod  l)  which  is  equivalent  to 


B'1  Nxn  =  B’b  (modi).  Substituting  in  this  constraint  and  relaxing  the  non¬ 
negativity  requirement  on  xB  yields 
Min  cNxN 

s.t.  B’'Nxn  =B''b  (modi)  (5.7) 

xN  >  0  and  integer 

Since  the  constraint  in  (5.7)  is  based  on  modulo  1,  all  that  is  needed  from  the 
columns  B'!N  and  B'lb  are  the  fractional  parts, =  B'1  and 

fi=B'1bj-\ji',b; J.  For  example,  the  fractional  part  of  3.25  is  0.25  and  the 

fractional  part  of  -3.25  is  0.75.  Let  ctj  be  a  column  vector  containing  the  fractional 
parts  of  B~'Nj  and  «0  be  a  column  vector  containing  the  fractional  parts  of  B'b.  The 
resulting  GMP  from  (5.1)  is 


Min  cNxN 

s-t.  ’Yj"j=\ajxj  =  ao  (mod  1) 
xN  >  0  and  integer 


(5.8  -  GMP) 


Example  5.1 

Consider  the  following  integer  linear  program  (ILP)  (Jensen  and  Bard  2003). 


Max  -2xj  +  x2 
s.t.  9xj-3x2>11 

X!  +  2x2<10  (PI) 

2x,  -  x2  <  7 
xpx2  >0  and  integer 
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The  feasible  region  is  shown  in  Figure  5.1.  The  optimal  solution  to  the  LP  relaxation 
is  X\  =  52/21,  x2  =  79/21,  and  z  =  -25/21.  The  dashed  lines  are  the  isovalue  lines  of  z 
and  the  shaded  area  is  the  feasible  region  for  PI  with  optimal  solution  x\  =  2,  x2=  2, 
and  z  =  -2,  which  is  not  the  closet  integer  point  to  x*p  in  Euclidean  distance.  It  is  the 
first  integer  point  encountered  as  we  move  the  isovalue  lines  of  z  into  the  feasible 
region. 


Adding  slacks  yields 

Max  -2x,  +  x2 

s.t.  9xj  -  3x2  —  Vj  =11 

xl+2x2  +s2  =  10 
2X[  -  x2  +  s3  =  7 

xx ,  x2 ,  s, ,  s2 ,  s3  >  0  and  integer 
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At  x’p  the  basic  variables  are  Vi,  jc2,  and  *3  with 


'9  -3  O' 

'  2/21 

3/21 

O' 

' 52/21 ' 

'-2/21 

3/21' 

B  = 

1  2  0 

B1  = 

-1/21 

9/21 

0 

B]b  = 

79/21 

B'N  = 

1/21 

9/21 

2  -1  1 

-5/21 

3/21 

1 

122/21] 

5/21 

3/21 

The  reduced  costs  for  *i  and  52  are  5/21  and  3/21  respectively.  The  resulting  GMP  is 

Min  5/21*,+ 3/2 1*2 

s.t.  19/21*!  +3/21*2  =10/21  (modi) 

1/21*! +9/21*2  =  16/21  (modi)  (GMP1) 

5/21*i  +3/21*2  =17/21  (modi) 

*1 ,  *2  >  0  and  integer 

Note  element  au  =  19/21  instead  of  -2/21  based  on  the  previous  definition  of 
fractional  part. 


5.1.2  The  Fractional  Group 

The  problem  defined  in  (5.8)  is  called  the  group  minimization  problem 
because  the  vectors  aj  and  «o  are  elements  of  a  finite  abelian  group  under  component¬ 
wise  addition  modulo  1.  The  group  equation  is  the  congruence  relationship 

Z”=i«7^=ao(modl)  (5‘9) 

Let  D  be  the  absolute  value  of  the  determinant  of  the  basis  B.  By  Cramer’s 
rule,  the  form  of  each  vector  element  ais  is  L  / D  where  Il}  is  integer  and  /,)  <  D. 

Since  each  element  of  the  vector  can  take  on  any  one  of  D  distinct  values,  the  size  of 

the  full  group  is  Dm  where  m  is  the  number  of  rows  in  the  problem  (size  of  the 

vector).  However,  the  elements  «,  generate  a  subgroup  of  size  <  D  (Johnson  1980, 
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Salkin  and  Mathur  1989).  Let  H(a)  =  {a ^  be  the  subgroup  generated  by  the 


columns  of  the  GMP. 
Example  5.2 

For  PI  we  have 


«!  =  Si  = 


19/21 

1/21 

5/21 


The  Iij  in  element  g  i  are  all  relatively  prime  with  D  =  21,  therefore  g i  generates  the  21 
element  cyclic  subgroup,  H(a )  (Salkin  and  Mathur  1989).  ( a2  generates  a  cyclic 
subgroup  of  order  7  which  is  subsumed  in  H(a).)  The  elements  of  H(a )  are  shown 


below. 


§1 

§2 

§3 

§4 

§5 

§6 

§7 

15/21 

3/21 

15/21^ 

RRI 

lnKC|l 

§8 

89 

§10 

§n 

§12 

§13 

§14 

■d£$K 

inS|| 

3/21 

9/21 

3/21 

"  1/21 " 
10/21 
8/21 

'20/2 f 
11/21 

1 3/  2 1  _ 

■ 

iKIAll 

§15 

§16 

§17 

§18 

§19 

§20 

§21 

12/21 

15/21 

12/21 

'  8/21 " 
17/21 
1/21 

'  6/21 ' 
18/21 
_6/21_ 

nlmal 

'  2/21 ' 
20/21 
16/21 

1 

1 
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We  have  a\  =  g  1,  and  «2  =  gf  =  gi*9  =  g9.  The  It]  in  element  g9  share  a 
common  divisor  with  D,  namely  3,  so  the  order  of  element  g9  is  21/3  =  7.  Since  ao  = 
g\\6  =  16gi  =  gi6,  setting  si  =  16  and  S2  =  0  satisfies  the  GMP1  constraints.  The 
corresponding  xb j  (from  Equation  5.4)  are  xi  =  4,  X2  =  3,  53  =  2.  This  yields  a  feasible 
but  not  optimal  solution  for  PI  with  z  =  -5.  A  lower  cost  solution  is  generated  by  oto 
=  gi6  =  gi  +  gis  =  gi  +  4gg  which  implies  that  sj  =  1  and  S2  =  4  and,  further,  that  xi  =  2, 
X2  =  2,  £3  =  5.  This  yields  the  optimal  solution  for  PI  with  z  =  -2  . 

Incidentally,  the  relaxations  of  the  rows  in  (5.8)  from  =  (modi)  to  >  are  the 

Gomory  fractional  cuts.  These  row  vectors  also  generate  a  finite  abelian  group  under 
component-wise  addition  modulo  1  that  has  a  size  of  <  D.  The  group  generated  by 
the  rows  and  the  group  generated  by  the  columns  are  isomorphic  (Gomory  1963). 
Gomory  also  notes  the  group  is  often  cyclic  and  therefore  can  be  generated  by  a 
single  element. 

5.1.3  The  Sufficient  Condition  for  xB  >  0 

The  GMP  attempts  to  find  the  integer  point  with  the  minimum  weighted  (by 
c  )  Euclidean  distance  to  x*p .  As  illustrated  in  Figure  5.1,  the  optimal  solution  to  an 

ILP  need  not  be  the  closest  integer  point  to  in  unweighted  Euclidean  distance. 
We  can  solve  the  GMP,  given  in  5.8,  for  any  ILP  and  when  xn  is  found,  xB*  can  be 
derived  using  Equation  5.4.  If  xB*  >0  then  x*  |xn*]  optimal  solution 
to  the  original  ILP. 
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Gomory  (1965,  1969)  provides  a  sufficient  condition  under  which  xb  will  be 
non-negative,  From  the  LP,  we  have  Bxb  =  b  and  xB  >  0 ,  which  implies  that  b  is 
contained  in  the  cone  generated  by  the  columns  of  B,  K{B).  Solving  the  GMP  implies 
an  associated  Xb-  Since  Bxb  =  b  -  Nxn,  if  b  -  Nxn  is  contained  in  K(B),  xB  >  0 . 

If  b  is  on  the  surface  of  K(B),  the  solution  is  degenerate  and  at  least  one  of  the 
basic  variables  is  equal  to  0.  Since  we  are  perturbing  h  by  (Nxn  )(. ,  we  need  to  have 
b  a  sufficient  distance  from  the  perimeter  of  K(B)  to  ensure  the  optimal  solution  to  the 
GMP  yields  xB  >0.  In  other  words,  the  basic  variable  values  from  the  LP  optimal 

must  be  sufficiently  greater  than  zero.  If  the  LP  yields  a  degenerate  solution,  then  b  is 
on  the  boundary  of  K(B).  Based  on  this  observation,  Balas  (1973)  shows  the 
sufficient  condition  is  never  satisfied  for  binary  programming  problems. 

5.1.4  Column  Reduction 

Gomory  (1969)  provides  two  cases  for  eliminating  a  column  from  the  GMP. 
From  a  strictly  GMP  perspective,  if  all  of  the  elements  in  a  non-basic  column  are 
integer,  that  redundant  column  will  map  to  the  group  identity  (zero  vector)  and 
contribute  nothing  to  forcing  the  basic  variables  to  be  integer.  Since  >0,  an  all 

integer  non-basic  column  will  never  appear  in  the  optimal  solution  to  the  GMP  and 
may  be  eliminated.  Similarly,  if  two  or  more  non-basic  columns  have  the  same 
fractional  components,  they  map  to  the  same  group  element.  Again,  from  a  strict 
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GMP  perspective,  the  column  with  the  smallest  c;.  can  be  retained  and  the  other 
redundant  columns  may  be  discarded. 

Unfortunately  these  redundant  columns  might  be  needed  to  construct  an 
optimal  solution  to  the  original  ILP.  If  the  optimal  GMP  solution  yields  a  negative 
j kb,  the  more  costly  “redundant”  columns  may  be  needed  in  the  ILP  to  ensure  xH  >  0 . 

For  example,  we  could  have  two  columns,  Nj  and  iV*,  with  identical  fractional 
components.  If  cy.  >  ck  and  ||al||  <  Jj/Vj ,  it  may  be  necessary  to  absorb  the  additional 
cost  of  c j  so  that  xB  >  0 .  Additionally,  if  we  have  xBi  <  0  and  an  all-integer  Nj  with 

<  0  we  may  need  to  use  Nj  to  ensure  xB  >  0 . 

Example  5.3 

Assume  the  solution  to  the  LP  relaxation  generated  the  following 


' 50/21 ' 

‘  21/21  8/21  29/21' 

cN  =  [2/21  3/21  2/21]  B'b  = 

110/21 

b'n  = 

42/21  5/21  26/21 

3/21 

-42/21  24/21  66/21 

The  associated  GMP  is 

Min  2/21xNl+3/2\xN2  +  2/2\xN3 
s.t.  0xNl+S/21xN2  +  S/2\xN3=S/2\  (modi) 

0xN1+5/2lxN2+5/2lxN3=5/2\  (modi) 
0xNl+3/2lxN1+3/2lxN3=3/2\  (modi) 
xNl,xN2,xN3  >  0  and  integer 
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The  first  column  has  all  integer  values  and  maps  to  the  group  identity.  The  second 
and  third  columns  have  the  identical  fractional  components  and  cN2  >  fv3 .  So  the 
first  and  second  columns  can  be  eliminated  from  the  GMP. 

Cleary  the  optimal  solution  to  the  GMP  is  xm  =  0,  xm  =  0,  and  =  1  with  a 
cost  of  2.  However,  this  GMP  solution  yields  xb\  -  1,  xgi  -  4,  and  xb3  =  -3  which  is 
an  infeasible  solution  to  ILP.  The  minimum  cost  solution  for  the  GMP  that  yields 
xB  >  0  is  xm  =  1,  xm  =  1,  and  xm  -  0  with  cost  5.  The  resulting  basic  variables  are 
xb\  =  1 ,  xb2  =  3,  and  xb3  =  1 . 


5.1.5  The  Factor  Group  Minimization  Problem  (FGMP) 

An  alternate  form  of  the  GMP  exists  based  on  a  factor  group  M(J)/M(B) 
where  M(I)  is  the  group  of  all  integer  vectors  of  size  m  and  M(B)  is  the  subgroup  of 
all  linear  integer  combinations  of  the  columns  in  the  optimal  LP  basis.  The  size  of 
the  factor  group  is  <  D  and  the  group  is  isomorphic  to  H(a),  therefore,  the  solutions  to 
both  versions  are  identical  (Gomory  1969,  Salkin  and  Mathur  1989). 

The  FGMP  formulation  is  found  by  finding  the  row  and  column  operations 
required  to  translate  B  to  the  Smith-Normal  form  and  performing  those  same  row  and 
column  operations  on  N  and  b.  If  the  method  used  to  solve  the  GMP  relies  on  the  full 
group  not  just  the  subgroup,  as  many  early  approaches  did  (Glover  1968,  Salkin  and 
Mathur  1989,  Wolsey  1971),  it  may  be  beneficial  to  use  the  FGMP  formulation 
(Johnson  1980). 
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5.1.6  Another  View  of  the  GMP 


We  can  view  the  GMP  generically  in  group  theoretic  terms.  Given  some 
arbitrary  group  G  we  want  to  find  the  minimum  cost  combination  of  group  elements 
that  sum  to  some  group  element  g0. 

Mi"  y.,,r/:,e>'lS> 

St.  IU«'<S>  =  So  <5-10> 

t(g)  >  0  and  integer  V  g  e  G 

where  c(g )  is  the  real-valued  weight  or  cost  for  group  element  g  and  t(g)  is  the  power 
or  multiplier  for  element  g  (Gomory  1969). 

5.2  Corner  Polyhedra 

In  this  section,  we  explore  the  geometry  behind  the  GMP.  Gomory  (1967, 
1969)  defines  corner  polyhedra  as  the  main  area  of  interest  in  an  ILP  and  the  feasible 
region  for  the  GMP. 

5.2.1  Corner  Polyhedra  in  X  Space 

Let  X  represent  the  original  decision  variables  of  the  ILP  (excluding  the  slack 
variables).  The  corner  polyhedron  in  X  space  (CP*)  is  found  by  first  relaxing  the 
bounds  on  the  basic  variables  from  the  optimal  solution  to  the  LP  relaxation  (Gomory 
1967,  1969).  For  each  Xj  e  X  ,  if  variable  Xj  is  basic  then  Xj  >  0  and  the  hyperplane 

Xj  =  0  does  not  pass  through  xLP" .  Similarly,  if  a  slack  variable  sj  is  basic  then  Sj  >  0 
and  age  <  £>„  where  a,  is  the  ith  row  of  A,  and  the  hyperplane  age  =  bt  also  does  not 
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pass  through  jclp*  .  So  relaxing  the  lower  bounds  on  the  basic  variables  is  equivalent 

to  relaxing  all  non-binding  constraints  at  *LP* . 

We  are  left  with  only  the  hyperplanes  passing  through  xLP’  which  form  an 

unbounded  polyhedral  cone  K.  All  of  the  feasible  integer  solutions  to  the  original  ILP 
are  contained  in  this  cone.  The  comer  polyhedron,  CP x,  is  the  convex  hull  of  the 
integer  points  in  K  (i.e.  CPX  =  conv(K)).  The  optimal  solution  to  the  ILP  is  a  vertex 
of  CP x  which  will  often  have  significantly  less  vertices  than  the  original  feasible 
region  conv(P)  (Gomory  1969). 

Example  5.4 

For  PI  at  xLP\  we  have  xh  x2,  and  S3  basic.  Relaxing  the  lower  bounds  on 

these  variables  yields  unbounded  polyhedral  cone  K  shown  with  solid  lines  in  Figure 
5.2.  The  associated  convex  hull  is  the  comer  polyhedron  shown  as  the  shaded  area  in 
Figure  5.2.  The  ten  original  feasible  points  are  labeled  A  through  J.  C  is  the  optimal 
ILP  solution. 
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*2  ▲ 


-  corner  polyhedron  in  x  space,  CPX 


Figure  5.2 


5.2.2  Corner  Polyhedra  in  XN  Space 

Let  Xn  represent  the  non-basic  variables  from  the  optimal  solution  to  the  LP 
relaxation,  the  variables  in  the  GMP  (including  non-basic  slack  variables).  We  can 

XN 

also  construct  the  corner  polyhedron  in  XN  space  ( CP  )  (Gomory  1969).  To  do  so 
we  relax  the  congruence  relationship  (5.9)  to 

Z"=1  auxj  -aio  Vi  =  i, ...  ro  (5.11) 

Certainly  any  feasible  point  for  (5.9)  satisfies  the  set  of  equations  in  (5.11).  The 

equations  in  (5.11)  are  the  classical  Gomory  fractional  cuts.  Let  P  be  the 
unbounded  convex  polyhedron  formed  by  (5.11)  and  the  non-negativity  of  x n-  The 

XN 

convex  hull  of  the  integer  points  in  P  that  satisfy  (5.9)  is  the  corner  polyhedron  in 
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XN  space,  CP  V  .  Any  vertex  v  of  CP x  corresponds  to  a  vertex  of  CP  N  and  vice- 


versa,  so  the  optimal  solution  to  the  original  ILP  is  also  a  vertex  of  CP  (Gomory 


Example  5.5 


we  have  si  and  S2  non-basic.  Relaxing  the  constraints  yields 


the  following  constraints 


The  convex  hull  satisfying  these  constraints  yields  the  corner  polyhedron  in  XN  space 


as  shown  as  the  shaded  area  in  Figure  5.3 


Figure  5.3  -  comer  polyhedron  in  Xn  space,  CP 
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The  reduced  cost  of  each  non-basic  variable  represents  the  cost  of  moving  one 
unit  away  from  the  origin  in  Xn  space  along  that  variable’s  axis.  Of  course,  not  every 
integer  point  in  Xn  space  yields  integer  basic  variables  when  xB  is  calculated  using 
Equation  5.4,  only  those  points  satisfying  (5.9).  These  integer-basis  points  are  circled 
in  Figure  5.3  and  labeled  based  on  their  corresponding  point  in  X  space.  For  instance, 
point  K  is  an  integer-basis  point  and  is  GMP1  feasible  but  it  is  infeasible  for  PI  as  it 
yields  xB  <  0.  The  integer-basis  points  comprise  a  sub-lattice  on  the  integer  lattice  of 

XN 

CP  .  Of  course  this  sub-lattice  is  not  random;  it  follows  a  distinct  pattern  based 

X„ 

on  the  group  H(a).  The  density  of  these  points  within  the  integer  lattice  of  CP  is 
related  to  the  order  of  the  group  elements  associated  with  each  non-basic  variable. 
Theorem  5.1 

The  frequency  of  integer-basis  points  along  axis  k  of  the  integer  lattice  is  order(a*). 
Proof 

Let  xN  be  an  integer-basis  point  in  Xy  space.  Fix  Xj  =x .,  Vxj  e  XN  then  choose 
some  xk  e  X N  to  become  free.  So  we  have  one  degree  of  freedom  and  are  exploring 
integer  lattice  points  along  k']l  axis.  From  the  congruence  relationship  (5.9)  we  know 
the  point  is  an  integer-basis  point  iff 

akxk  =ao~ ajxj  =  a  (mod  l) . 

j*k 

We  know  xk  =  xk  satisfies  the  above  relationship  and  therefore  a  is  a  power  of  a*.  If 
akxk  =a'(mod  l),  then  we  must  have  ak  (x*  +  p  order  (oq))  =  a  (mod  1)  for  all 
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integer  p  by  the  definition  of  the  order  of  a  group  element.  Therefore  the  frequency 
of  integer-basis  points  along  axis  k  of  the  integer  lattice  is  order(a^).  q> 

As  noted  earlier,  the  form  of  the  elements  of  o.j  is  / . . / D  where  ly  is  integer 

and  Iy  <  D.  Let  gcd (D,/7)  be  the  greatest  common  divisor  of  D  and  the  Iy  over  j.  The 
order  of  element  «,  is  D/gcd(D,Ij)  (Salkin  and  Mathur  1989). 

Example  5.6 

In  problem  PI,  we  have  order(«i)  =  21  and  order^)  =  7.  In  Figure  5.3,  we 
can  clearly  see  that  integer-basis  points  are  21  lattice  points  apart  along  the  s\  axis 
and  7  lattice  points  along  the  S2  axis. 

The  order  of  any  element  in  a  group  must  be  less  than  the  size  of  the  group. 
Since  the  size  of  H{a)  is  <D,  D  gives  us  an  upper  bound  on  the  frequency  of  integer- 
basis  points  along  any  axis.  It  is  well  known  that  if  D  -  1,  the  LP  relaxation  will 
yield  an  all  integer  basis.  In  GMP  terms,  D  -  1  implies  || H  («)|  =  1  and  every  integer 

point  in  Xn  space  is  an  integer-basis  point. 

If  B]b  is  integer  then  «o  in  (5.8)  is  the  group  identity  or  zero  vector.  In  this 
case,  xn  =  0  is  a  feasible  solution  to  the  GMP  and  the  origin  in  Xn  space  is  on  the 
integer-basis  sub-lattice  and  inside  the  corner  polyhedron.  Since  the  group  identity 
must  be  a  power  of  every  group  element,  there  are  integer-basis  sub-lattice  points  on 
each  axis  in  Xn  space.  So  the  comer  polyhedron  becomes  the  positive  quadrant  and 
the  origin  its  only  vertex. 
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5.3  Local  Search  Neighborhoods  in  XN  Space 

We  can  solve  the  GMP  by  searching  the  corner  polyhedron  in  XN  space.  One 
simple  neighborhood  is  to  increment  or  decrement  by  1  a  single  non-basic  variable 
xnj  with  cost  ±Cj  for  j  =  1  to  n.  This  neighborhood  would  traverse  the  integer  lattice 

in  Xn  but  would  require  a  penalty  function  to  drive  the  search  to  the  sub-lattice  points 
yielding  an  integer  basis.  While  this  neighborhood  is  simple  to  implement  and  allows 
us  to  take  advantage  of  the  fine  grain  gradient  provided  by  the  individual  c. ,  it  is 

affected  by  the  size  of  D.  Experimentation  has  shown  it  to  be  ineffective  for  all  but 
very  small  problems.  The  distance  between  the  points  on  the  sub-lattice  is  too  great 
and  as  the  problems  get  larger  the  algorithm  failed  to  find  a  single  sub-lattice  point. 
So  let  us  consider  how  to  traverse  the  integer-basis  sub-lattice. 

5.3.1  Identity  Moves 

Let  xN  and  yN  be  integer-basis  points.  Since  both  satisfy  the  congruence 

relationship  (5.9),  define  d  =  xN  -  yN ,  so  Z  ■=!«/,  =0(mOdl)  with  0  being  the 

group  identity.  Note  that  xN  +  p  d  is  an  integer-basis  point  for  all  integer  p.  The  n- 
dimensional  vector  d  allows  us  to  move  from  one  integer-basis  point  to  another  in 
direction  of  d;  it  captures  the  direction  and  the  step  size. 

The  move  value  cNd  measures  the  change  in  weighted  Euclidean  distance  to 

xLP’  and  is  independent  of  the  current  point.  Given  a  known  set  of  identity  moves, 
we  can  order  them  by  move  value  from  best  to  worst.  At  every  iteration  we  take  the 
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first  available  move  based  on  our  search  parameters  and  avoid  searching  the  entire 
neighborhood. 

Example  5.7 

In  Figure  5.4,  =  (13,  5)  with  zgmp  =  80/21  and  J n  =  (16,  7)  with  zgmp  = 


101/21.  Letting  d  =  Jn  -  Gn  =  (3,  2),  then  cNd  =  21/21.  Moving  along  this  line  we 

have  G/v  -  d  =  (10,  3)  =  DiV  with  zgmp  =  59/21  and  Gjv  -  2 d=  (7,  1)  =  An  with  zgmp  = 
38/21.  Calculating  xg  using  Equation  5.4  yields  solutions  D  and  A  from  Figure  5.2 
with  zip  =  -4  and  zip  =  -3  respectively. 


1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  s i 


Figure  5.4  -  identity  moves  in  Xn  space 


Deriving  identity  moves  as  the  difference  between  two  solutions  does  not 
require  either  solution  to  be  feasible.  The  solutions  of  course  must  be  integer. 
However,  they  may  have  some  Xb  <  0  and/or  xn  <  0.  Since  we  calculate  Xg  from 
Equation  5.4  will  have  always  have  Ax  =  b  .  The  difference  between  two  infeasible 
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moves  still  gives  us  a  valid  identity  move  and  we  may  be  able  to  use  such  moves  to 
guide  the  search  into  the  feasible  region. 

Example  5.8 

Although  solutions  x  =  (0,  0,  -11,  10,  7)f  and  y  =  (0,  1,  -14,  8,  8/  satisfy 
Ax  =  b  from  PI,  they  are  infeasible  because  both  have  Sj  <  0.  The  associated  non- 
basic  variable  vectors  arejc;v  =  (-11,  10)  andy;v  =  (-14,  8).  Letting  d  =  : =  (3,  2) 
gives  us  the  same  identity  move  found  earlier.  Unfortunately,  since  the  vector 
connecting  x\  and  y,v  does  not  intersect  the  PI  feasible  region,  we  cannot  reach  the 
feasible  region  using  just  this  move. 

5.3.2  Generating  Identity  Moves 

Generating  identity  moves  as  the  difference  between  integer  solutions  can  be 
difficult  if  integer  solutions  are  hard  to  come  by.  Fortunately,  these  solutions  need 
not  be  feasible.  We  now  define  a  set  of  easily  generated  atomic  identity  moves.  As 
Theorem  5.3  will  show,  all  possible  identity  moves  can  be  expressed  as  a  linear 
integer  combination  of  these  atomic  moves. 

Using  Xj  =  0  e  X  as  our  baseline  solution,  we  can  create  additional 

solutions  by  setting  each  Xj  to  1  in  turn  and  adjusting  the  slack  variables.  We  then 
subtract  the  difference  between  each  new  solution  and  the  baseline  solution 
(including  corresponding  slack  variables).  The  resulting  difference  vector  details  the 
change  in  each  non-basic  variable  and  each  basic  variable  as  a  result  of  the  move. 
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The  value  of  the  move  can  be  calculated  by  multiplying  the  vector  of  reduced  costs 


by  the  difference  vector  since  c }  -  0  if  xj  is  basic. 

Here  is  the  procedure  that  builds  atomic  identity  move  j  (denoted  d1)  for  j  <  n: 

1.  Set  Vxj  =  1,  Vxk  =  0  \fxk  eX,k*  j 

2.  For  each  slack  variable  ^ 

f -ar  if  constraint  i  is  < 

Vs,  =  1  .  ,y  .  ..  where  an  is  from  A  in  Equation  5.1 

Ici'j  if  constraint  i  is  >  1 

3.  d1  =  [Vx|VsJ  has  cost  cdJ 
Example  5.9 

For  PI,  build  rf1  from  Vx{  =  1 ,  Vx2  =  0 ,  Vs,  =  9  ,  Vs2  =  -1 ,  and  Vs3  =  -2  so 

d1  =  (1,  0,  9,  -1,  -2)f.  The  reduced  cost  vector  is  c  =  (0,  0,  5/21,  3/21,  0)  so  erf1  =  2. 
For  rf2,  use  Vx,  =  0 ,  Vx2  =  1 ,  Vsj  =  -3 ,  Vs2  =  -2 ,  and  Vs3  =1  so  rf2  =  (0,  1 ,  -3,  -2,  1)' 

and  erf2  =  -1.  Starting  from  any  integer-basis  point  in  XN  space  we  can  add  linear- 
integer  combinations  of  these  moves  and  reach  another  integer-basis  point.  Going 
back  to  our  previous  example  x  =  (0,  0,  -11,  10,  7)(  with  ex  =  0.  We  add  the  two 
identity  moves  to  get  x  +  2rf'  +  2 rf2  =  (2,  2,  1,4,  5/  =  C  and  ex  -  2crf1  -  2 erf2  =0-4 
+  2  =  -2. 

Since  the  rf’ s  change  basic  and  non-basic  variables  to  guarantee  an  integer 
basis,  we  need  not  compute  the  basic  variables  using  Equation  5.4.  We  can  simply 
update  Xb  using  the  rf’s.  The  atomic  identity  move  d1  is  column  j  of  the  ( n+m )  by  n 

matrix  whose  upper  part  is  the  n  by  n  identity  and  whose  lower  part  is  VA  where  V  is 
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the  diagonal  matrix  of  constraint  signs  (1  for  >  and  -1  for  <)  and  A  is  from  Equation 
5.1. 

Theorem  5.2 

The  cost  of  d 1  is  Cj . 

Proof 

The  cost  of  d1  is 

cdJ  =E-[Vx|VxJ  =cj+Yifn+i  -V.v,  (5.13) 

If  constraint  i  is  >  then  cn+i  =  nt  and  Vs,  =  atj .  If  constraint  i  is  <  then  cn+:  =  -ni  and 
Vs,  =  -atj .  So  (5.13)  becomes 

cdJ  =c.  +Ym*,  -a..  =  c-  +jiA =c,  7i A .  +nA.  =  c.  □ 

J  1  1  U  J  J  J  J  J  J 

Theorem  5.2  should  not  be  surprising  as  dj  only  changes  one  variable  with  a  non-zero 
cost,  xj.  So  when  we  traverse  the  integer-basis  sub-lattice  in  Xn  space  via  identity 
moves,  we  are  simultaneously  exploring  the  integer  lattice  in  X  space  and  vice-versa. 
By  viewing  the  search  in  terms  of  XN  space,  we  gain  insight  we  can  use  to  enhance 
our  search. 

Example  5.10 

Let  Gjv  =  (13,  5)  be  the  incumbent  solution.  Using  Equation  5.4  yields  xb  =  (3,  1,  2) 
so  Gn  is  the  Xn  space  image  of  solution  G  in  X  space  (Figure  5.2).  The  cost  of 
solutions  G  and  Gn  is  -5.  If  we  apply  move  d2  =  (0,  1,  -3,  -2,  l)f,  we  move  from  Gn  = 
(13,  5)  to  Dn  =  (10,  3).  Using  Equation  5.4  yields  xB  -  (3,  2,  3)  so  Dn  is  the  Xn 
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space  image  of  solution  D  in  X  space  (Figure  5.2).  The  cost  of  solutions  D  and  Dn  is 
-4.  Move  d 2  increased  xt  by  1  and  did  not  change  x\  so  the  change  in  solution  value 
was  C2  =  1  as  expected. 

Theorem  5.3 

Any  move  from  one  integer-basis  point  to  another  can  be  expressed  as  an  integer 
linear  combination  of  the  atomic  identity  moves. 

Proof 

There  is  a  bijection  between  the  integer  points  in  X  space  and  the  integer-basis  points 
in  Xn  space.  Given  any  integer  point  in  X  space,  we  can  calculate  the  corresponding 
slack  variables.  Removing  the  basic  variables  gives  us  the  corresponding  point  in  X n 
space.  Given  a  point  in  XN  space,  we  can  calculate  the  basic  variables  using  Equation 
5.4.  Extracting  the  variables  in  X  gives  us  the  corresponding  point  in  X  space.  Given 
any  two  points  in  X  space,  the  difference  between  them  is  an  n-dimensional  vector 
that  can  be  expressed  as  an  integer  linear  combination  of  n  rc-dimensional  unit 
vectors.  We  have  n  atomic  identity  moves  each  changing  a  single  xk  e  X  by  1  while 

fixing  all  other  Xj  el  at  their  current  values  so  move  d 1  is  equivalent  to  the  n- 

dimensional  unit  vector  j  in  X  space  and  the  same  integer  linear  combination  of 
atomic  identity  moves  transitions  between  the  images  of  the  two  points  in  Xn  space.  □ 

Also  by  Theorem  5.3,  the  atomic  identity  moves  provide  us  with  connectivity  in  our 
search  space. 
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5.3.3  Identity  Move  Neighborhoods 

The  search  neighborhood  is  comprised  all  of  the  solutions  reachable  from  the 
current  solution  in  one  iteration.  Given  a  current  integer  solution  xr,  a  move  from  the 
current  solution  to  the  next  is  defined  as 

*'+l=*r+2"=i Pjd'  Pj£  V/  =  l...n  (5.14) 

However,  evaluating  all  such  combinations  is  equivalent  to  enumeration.  A 
larger  neighborhood  increases  the  likelihood  of  finding  good  solutions  but  also 
increases  the  computational  effort  required,  so  we  must  compromise.  We  can  restrict 
the  neighborhood  size  by  imposing  a  limit  on  the  atomic  identity  move  coefficients. 

Pi- 

We  define  a  k-step  move  as  one  where  =  then  k  restricts  our 

neighborhood  size.  We  can  have  multiple  levels  of  k  where  &i+&2-step  implies  the 
restriction  =&,  or  =  k2 .  We  can  further  restrict  the  neighborhood 

size  by  imposing  bounds  on  the  individual  coefficients  as  l<Pj<  u  .  Of  course  if 

\l\>k  and/or  \u\>k  then  k  becomes  the  only  restriction.  Similarly  if  £>«|/|  and 

&>nj«|  then  /  and  u  provide  the  only  restrictions.  Finally,  we  can  denote  an 

unrestricted  parameter  as  the  wildcard  *.  Therefore  a  neighborhood  N  can  be  defined 
by  (5.14)  and  the  triplet  ( k ,  l,  u ).  For  example,  the  neighborhood  iV(l,-l,l)  is 
comprised  of  all  solutions  reached  by  adding  or  subtracting  a  single  atomic  identity 
move,  the  1-step  identity  move  neighborhood.  Regardless  of  the  neighborhood 
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chosen,  the  costs  of  the  moves  are  still  fixed.  Therefore  we  can  sort  them  from  best 
to  worst  and  choose  the  first  feasible  move. 

5.3.4  Bound  Strengthening 

We  can  use  the  reduced  costs  and  the  geometry  of  XN  space  to  strengthen  the 
bounds  on  the  non-basic  variables  as  the  search  progresses.  Let  |GMP  be  the  objective 

value  of  the  current  best  solution  found  to  the  GMP,  jrGMP .  Then  £GMP  is  the  gap 
between  the  objective  value  of  the  corresponding  integer  solution  to  the  ILP  and  the 
objective  value  to  the  optimal  solution  of  the  LP  relaxation,  zu* .  If  c .  >  zGMP  then 
non-basic  variable  Xj  must  be  0  in  any  solution  better  than  iGMP .  Otherwise  Xj  can  be 
nonzero  in  solutions  better  than  where  [loMp/cj.J  bounds  Xj. 

5.4  Conclusions 

The  combinations  of  non-basic  variables  that  yield  an  integer  basic  variables 
form  a  sub-lattice  in  X #  space.  The  density  of  this  sub-lattice  is  bounded  by  the 
determinant  of  the  LP  basis.  Any  search  of  the  integer  feasible  sub-lattice  in  Xn 
space  must  incorporate  the  atomic  identity  moves.  Using  the  full  set  of  moves 
provides  us  with  connectivity  throughout  the  solution  space  (Theorem  5.3).  The 
neighborhoods  presented  here  are  general  enough  to  be  applied  to  almost  any 
metaheuristic  search  method.  We  present  an  example  Group  Theoretic  Tabu  Search 
implementation  of  these  ideas  in  the  next  chapter. 
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Chapter  6  -  A  Group  Theoretic  Tabu  Search  Algorithm  for 
the  Group  Minimization  Problem  (GTTS-GMP) 

In  this  chapter  we  develop  a  group  theoretic  tabu  search  (GTTS)  algorithm  for 
solving  the  group  minimization  problem  (GMP),  the  GTTS-GMP.  Given  an  integer 
linear  program  (ILP),  we  solve  a  linear  programming  (LP)  relaxation  of  the  ILP  and 
use  the  LP  optimum  solution  to  formulate  the  GMP.  We  use  identity  move 
neighborhoods  to  explore  the  solution  space  and  bound  strengthening  based  on 
reduced  costs  as  described  in  Section  5.3.4. 

As  previously  noted,  an  algorithm  that  solves  the  GMP  solves  the  general  ILP. 
Tabu  search  implementations  for  the  general  ILP  have  been  largely  ignored. 
Algorithms  developed  to  exploit  the  specific  characteristics  of  special  case  problems 
tend  to  be  more  effective  and  dominate  the  literature.  However,  a  general  purpose 
implementation  is  important  to  develop  techniques  that  are  universal  in  their 
effectiveness  and  to  avoid  developing  a  new  implementation  for  each  new  problem 
set  (Glover  and  Laguna  1997). 

We  use  multi-dimensional  knapsack  problems  (MDKP),  both  integer  and 
binary,  and  set  covering  problems  (SCP),  both  unicost  and  weighted,  to  test  the 
versatility  of  our  approach.  GTTS-GMP  performs  very  well  on  this  diverse  problem 
set. 
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6.1  Methodology 

We  will  first  present  a  short  overview  of  the  methodology  and  provide  details 
on  each  aspect  of  the  algorithm.  At  certain  points,  such  as  in  the  discussion  of  escape 
procedures,  more  than  one  alternative  will  be  described.  We  will  test  each  of  these 
alternatives  and  present  the  results  in  Section  6.2. 

6.1.1  Overview 

First,  CPLEX  solves  the  linear  relaxation  of  the  ILP.  The  associated  LP 
optimal  basic  variables,  xBi ,  and  non-basic  variables,  xNj ,  and  the  reduced  costs,  cN 

are  used  to  formulate  the  GMP.  A  GTTS-GMP  starting  solution,  x° ,  is  generated  by 
rounding  the  non-slack  xB  without  regard  to  constraint  feasibility.  The  start  values 
for  the  slack  variables,  both  basic  and  non-basic,  are  then  calculated  from  Ax  and  b. 
A  l-step  identity  move  and  its  inverse  is  created  for  each  of  the  original  decision 
variables  and  placed  in  move  list  in  ascending  order  by  cost.  A  reactive  tabu  search 
(RTS)  procedure  is  used  to  explore  the  solution  space  using  an  identity  move 
neighborhood.  When  cycling  is  detected  an  escape  procedure  is  invoked.  The 
algorithm  terminates  when  a  GTTS-GMP  solution  value  equal  to  min(c  )  is  found  or 

when  the  allotted  time  has  expired. 
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6.1.2  Solving  the  LP  Relaxation 

The  pseudo  code  for  the  initialization  phase  of  the  algorithm  is  shown  in 
Figure  6.1.  CPLEX  is  initialized  and  the  problem  data  is  input.  We  extract  the 
problem  information  from  CPLEX  and  relax  the  integrality  restriction  on  the  decision 
variables.  CPLEX  is  then  asked  to  solve  the  problem  and  the  solution  and  reduced 
cost  vectors  are  retrieved. 

In  the  previous  chapter,  we  assumed  the  original  ILP  was  a  maximization 
problem  with  no  upper  bounded  variables.  When  we  remove  these  assumptions,  we 
introduce  some  reduced  costs  that  are  not  positive  and  some  non-basic  variables  that 
are  not  zero  (at  their  lower  bound)  in  the  LP  optimal  solution.  If  the  original  problem 
is  a  maximization  problem,  a  non-basic  variable  at  its  upper  bound  will  have  cj  <  0 

while  a  non-basic  variable  at  its  lower  bound  will  have  c .  >  0  .  If  the  variable  is  at  its 
upper  bound,  cj  represents  the  penalty  for  decreasing  the  variable  by  1.  If  the 
variable  is  at  its  lower  bound,  ~c,  represents  the  penalty  for  increasing  the  variable  by 
1.  Let  V*.  be  the  change  in  non-basic  variable  j.  Our  objective  function  value 
zgmp  =  ^Cj  V*y.  regardless  of  which  variables  are  at  their  upper  or  lower  bounds. 
For  a  minimization  problem,  the  signs  on  cj  are  reversed,  but  zgmp  is  calculated  in 
the  same  manner.  For  maximization  zGMP  >  0  and  for  minimization  zGMP  <  0 .  In 
both  cases,  Zjp  -^lp  ^gmp  * 
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The  reduced  costs  returned  by  CPLEX  are  prone  to  round-off  and  truncation 
error  and  occasionally  can  cause  a  small  valued  cj.  to  carry  the  wrong  sign.  For  this 

reason,  care  must  be  taken  to  avoid  exacerbating  this  error  when  using  the  cj .  As 

discussed  in  Section  5.3.2,  to  reduce  the  potential  for  error,  we  use  the  original 
objective  function  cost  c,  for  the  move  value  of  each  1-step  identity  move  (Theorem 
5.3). 


Relax  integrality  and  solve  the  LP  with  CPLEX 

Retrieve  the  reduced  cost  vector,  c  ,  optimal  LP  solution,  xLP* , 

optimal  LP  objective  function  value,  zjjp ,  and  basis  information 
Generate  the  starting  solution  jc° 

Round  each  non-slack  basic  variable  to  the  nearest  integer 
Calculate  the  slack  variables  as  s  =  b- Ax 
Calculate  the  objective  function 


if  max  z, 
else 


GMP 


=  zLp  ~cx 


,0 


^GMP  ^LP  +  CX 

Generate  1-step  identity  move  neighborhood,  N(  1,  -1,  1) 


Figure  6.1  -  Initialization  Phase 


6.1.3  Finding  a  Starting  Solution 

To  begin  our  search,  we  must  have  a  starting  solution.  Generating  a  feasible 
starting  solution  for  a  general  IP  is  not  an  easy  task.  Typically,  problem  specifics 
must  be  used  to  create  a  feasible  starting  solution  efficiently.  We  have  shown  in  the 
previous  chapter  the  use  of  identity  moves  allows  us  to  move  from  a  solution  outside 
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the  feasible  region  to  a  solution  inside  the  feasible  region.  So,  for  our  search 
methodology,  feasibility  of  the  starting  solution  is  not  as  important  as  integrality. 

We  would  also  like  to  be  close  to  xLP’  to  hopefully  reduce  the  number  of 
iterations  required  to  move  to  an  area  containing  good  solutions.  Perhaps  the  easiest 
way  to  do  this  is  to  round  xLP* .  We  round  each  non-slack  basic  variable  to  the 


nearest  integer  and  then  recalculate  the  slack  variables  as  s  =  b-  Ax  .  This  produces 


an  integer  solution  on  a  vertex  of  the  hypercube  surrounding  xLP*  in  X  space.  The 
solution  may  or  may  not  be  feasible. 


6.1.4  Identity  Move  neighborhoods 

Theorem  5.3  assures  that  any  neighborhood  containing  the  1-step  identity 
move  neighborhood,  iV(l,-l,l),  will  guarantee  connectivity  throughout  the  solution 
space.  Futher,  Theorem  5.3  assures  any  identity  move  neighborhood  is  a  composite 
of  iV(l,-l,l).  However,  larger  composite  move  neighborhoods  may  perform  better 
because  they  reach  more  solutions  in  a  single  iteration.  The  size  of  iV(l,-l,l)  is 
simply  2 n. 

The  full  2-step  identity  move  neighborhood,  N( 2, -2, 2),  contains  all  solutions 
reached  from  the  current  by  combining  any  2  1-step  moves.  The  size  of  N( 2, -2,2)  is 


Since  a  small  problem  with  n 


1000  implies  a 


neighborhood  with  precisely  2  million  members,  using  N(2,-2,2)  is  not  practical  for 
most  problems. 


Candidate  lists  are  often  used  in  tabu  search  to  restrict  neighborhood  size 
(Glover  and  Laguna  1997).  For  example,  in  a  TSP  we  may  restrict  swap  moves  to 
only  consider  swapping  cities  within  ten  positions  of  each  other  in  the  current 
solution.  Such  a  restriction  is  often  arbitrary,  but  is  valid  as  long  as  connectivity 
within  the  solution  space  is  maintained.  Using  Ml, -1,1)  in  conjunction  with  a 
candidate  list  comprised  of  any  subset  of  N( 2,-2, 2)  assures  connectivity. 

One  possible  candidate  list  is  constructed  as  follows:  create  JV(1,-1,1)  and  sort 
the  moves  in  ascending  order  by  move  value.  (This  list  is  symmetric:  the  bottom  of 
the  list  is  the  inverse  of  the  top.)  Next  combine  each  move  with  the  2  moves 
immediately  preceding  and  following  its  inverse  in  the  list.  This  strategy  excludes  the 
extreme  combinations,  combining  two  really  good  moves  or  two  really  bad  moves, 
and  creates  “fine-tuning”  moves  with  smaller  move  values.  The  number  of  combined 
moves  is  2 n  -  2,  so  the  size  of  the  neighborhood  is  2n  +  2n  -  2  =  An  -  2. 

Example  6.1 


Ml, -1,1) 


Move  3 

-150 

Move  5 

-140 

Move  2 

-115 

Move  1 

-90 

Move  4 

-80 

Move  -4 

80 

Move  -1 

90 

Move  -2 

115 

Move  -5 

140 

Move  -3 

150 

Subset  of  M2, -2, 2) 


Move  5,  -2 

-25 

Move  2,-1 

-25 

Move  3,  -5 

-10 

Move  1 ,  -4 

-10 

Move  4,-1 

10 

Move  5,  -3 

10 

Move  1,  -2 

25 

Move  2,  -5 

25 
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An  alternate  approach  is  to  create  combination  moves  dynamically  during  the 
search.  Given  that  the  best  feasible  non-tabu  move  is  taken  at  each  iteration  and  that 
the  values  of  the  moves  are  fixed,  if  the  move  just  executed  is  feasible  it  will  be 
chosen  again.  This  will  continue  until  the  move  is  no  longer  feasible.  In  other  words, 
we  have  reached  a  boundary  of  the  feasible  region.  So  changes  in  a  chosen  move 
occur  at  the  boundaries  of  the  feasible  region. 

We  can  create  iV(l,-l,l)  as  described  above  and  begin  our  search.  When  the 
current  1-step  move  being  executed  is  different  than  the  previous  one,  we  combine 
the  two  moves  and  add  the  new  composite  move  (and  its  inverse)  to  the  move  list. 
Our  algorithm  now  learns  the  “shape”  of  the  solution  space.  Of  course,  a  restriction 
will  need  to  be  in  place  to  limit  the  number  of  dynamic  moves  created  in  this  manner. 

Regardless  of  the  neighborhood  used,  the  moves  are  sorted  in  ascending  order 
and  the  first  (best)  feasible  non-tabu  move  is  chosen.  If  a  feasible  non-tabu  move 
cannot  be  found,  the  first  infeasible  non-tabu  move  is  chosen.  If  the  current  solution 
is  infeasible,  the  first  non-tabu  move  reducing  infeasibility  is  chosen.  If  a  non-tabu 
move  reducing  feasibility  cannot  be  found,  the  first  non-tabu  move  is  chosen. 

Finally,  if  two  or  move  moves  have  the  same  value  and  are  non-tabu  the  move 
with  the  greatest  normalized-surplus  is  chosen.  The  normalized  surplus  is  sum  of  the 
surplus  in  the  constraints,  each  normalized  by  their  respective  right  hand  side  values. 
We  are  essentially  choosing  the  more  “interior”  solution.  Any  further  ties  are  broken 
lexicographically.  In  our  testing,  we  compare  the  dynamic  approach  versus  the  basic 
A(l,-l,l)  neighborhood. 
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6.1.5  Tabu  Search  and  Tabu  Structures 

The  pseudo  code  for  the  search  phase  is  shown  in  Figure  6.2.  When  using  a 
search  neighborhood,  we  may  become  trapped  in  local  optima.  Tabu  search  (TS) 
helps  us  to  escape  local  optima  by  allowing  non-improving  moves  and  making 
recently  visited  solutions  tabu  (Glover  and  Laguna  1997). 

After  a  move  is  executed  an  attribute  of  the  current  solution  or  move  is  added 
to  the  tabu  list.  For  a  number  of  iterations,  tabu  tenure,  solutions  or  moves  on  that 
list  are  not  permitted  unless  such  a  move  leads  to  a  solution  superior  to  any  found  to 
that  point  in  the  search.  Solution  attribute  strategies  typically  allow  a  more  flexible 
search  than  move  based  strategies. 

From  Section  5.3.3,  a  move  neighborhood  is  specified  by  the  triplet  ( k ,  l,  u) 
where  k  is  the  sum  of  the  multipliers  for  the  atomic  moves,  /  is  the  lower  bound  of  the 
multipliers  and  u  is  the  upper  bound  of  the  multipliers.  If  l  =  -u  then  we  have  a 
symmetric  neighborhood.  When  a  move  is  performed  we  can  add  its  inverse  move  to 
the  list.  This  tabu  strategy  would  work  for  simple  neighborhoods  like  N(  1,  -1,  1),  but 
for  more  complex  neighborhoods  where  k  >  1  there  will  be  more  than  one  way  to 
return  back  to  the  same  solution.  For  such  neighborhoods  a  solution  attribute  is 
required  if  returning  to  a  previously  visited  solution  within  its  tabu  tenure  is  to  be 
prohibited. 

In  a  reactive  tabu  search  (RTS)  we  also  implement  a  tabu  structure  to  track  the 
frequency  in  which  solutions  are  visited  (Battiti  and  Tecchiolli  1994).  This  structure 
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is  used  to  detect  the  whether  the  search  is  trapped  in  an  attractor  basin  and  to  adjust 
the  tabu  tenure  or  implement  escape  procedures  when  such  a  basin  is  detected.  When 
a  solution  is  repeated  within  a  predetermined  interval  the  tenure  is  increased  to 
facilitate  diversification  into  new  areas  of  the  solution  space.  When  a  predetermined 
number  of  consecutive  new  solutions  are  encountered,  tenure  is  reduced  to  facilitate 
intensification  of  the  search  within  the  current  subset  of  the  solution  space.  When  a 
predetermined  number  of  solutions  have  been  repeated,  we  conclude  the  presence  of 
an  attractor  basin  and  an  escape  procedure  is  invoked. 

Due  to  the  immense  size  of  the  solution  space,  we  can  not  explicitly  store  each 
solution  encountered.  Further,  since  all  previous  solutions  must  be  checked  at  every 
iteration  we  must  have  an  efficient  way  to  store  and  access  solutions.  Just  as  in 
Section  4.2.5,  a  hash  function  (Woodruff  and  Zemel  1993),  <p(x) ,  is  used.  A  random 
number,  Pj,  is  generated  for  each  x The  solution’s  hash  value  is 

<P{x)  =  l[j‘j-\Pjxj  •  Solution  information,  such  as  number  of  visits  and  iteration  last 

visited,  is  then  stored  in  a  table  based  on  (p{x) .  When  two  or  more  distinct  solutions 

produce  the  same  hash  value,  a  collision  occurs.  We  can  limit  collisions  by  further 
distinguishing  a  solution  by  its  deficit  or  surplus. 

We  can  implement  a  solution-based  tabu  strategy  using  the  solution  tabu 
structure  and  <p(x).  Since  (p(x)  is  an  additive  function,  we  can  calculate  the  change 

in  #?(x),  A^?(x),  for  each  move.  As  we  evaluate  each  move,  we  can  use  A^(jc), 
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the  solution  tabu  structure  and  the  current  solutions  hash  value  to  determine  if  the 
move  leads  to  a  recently  visited  solution. 

We  use  a  combination  of  move-based  and  solution-based  tabu  strategies  in  our 
algorithm.  A  move  is  tabu  if  the  move’s  inverse  has  been  executed  in 
MOVE_TENURE  iterations  or  if  the  solution  has  been  visited  in  SOL_TENURE 
iterations.  If  a  move  leads  to  a  new  best  found  solution,  tabu  status  is  ignored. 


Search  loop  iteration  i 
If  jr'  is  feasible 

Find  the  first  feasible  non-tabu  move 

else 

Find  the  first  non-tabu  move  the  reduces  infeasibility 
Execute  the  move 

If  jc'+1  feasible  and  better  than  best  found 
Update  best  found 
Check  for  optimality 
Update  bounds 

If  using  dynamic  neighborhoods 

If  the  current  move  is  different  then  the  previous  one 

Combine  the  moves  and  add  the  new  move  to  the  list 
Update  move  and  solution  tabu  structures  and  get  search  status 
If  status  is  intensification 

Decrease  tabu  tenure 
If  status  is  diversification 
Increase  tabu  tenure 
If  status  is  cycling 
Escape 

If  stopping  criteria  not  met  repeat  search  loop 
Return  best  solution  found 


Figure  6.2  -  Search  Phase 
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6.1.6  Escape  Procedures 

We  considered  two  escape  strategies  -  random  restart  and  using  a  surrogate 
objective  function.  Both  escape  strategies  strive  to  break  the  current  cycle  or  capture 
basin  while  returning  the  search  to  the  vicinity  of  xLP* . 

For  random  restart,  we  round  xLP*  in  a  random  fashion.  For  each  non-integer 

Xj  the  probability  it  is  rounded  to  |~ Xj  "j  is  |~jty. "]  -  x}  and  the  probability  it  is  rounded  to 

is  The  current  solution  is  then  set  to  the  newly  rounded  xLP* 

solution.  Again  this  solution  may  or  may  not  be  feasible.  The  move-based  tabu 
information  is  reset  as  it  was  based  on  the  previous  location  and  the  search  is 
resumed.  When  generating  a  random  solution,  we  do  not  allow  a  restart  solution  to 
be  a  solution  that  has  been  visited  more  than  MAX_RESTART_REPEAT  times. 

For  the  surrogate  objective  function  strategy,  we  change  our  objective  to  that 
of  minimizing  the  Euclidean  distance  to  xLP* .  The  surrogate  objective  function  is 

i;,h-4v)2  (6-1) 

Minimizing  Equation  6.1  requires  us  to  search  the  entire  neighborhood  not  just  the 
top  of  the  sorted  list.  We  do  require  the  Xj  to  be  within  their  bounds  but  do  not  require 
the  solution  to  be  feasible  with  respect  to  the  constraints. 

Of  course,  an  optimal  solution  in  terms  of  Equation  6.1  is  our  starting 
solution,  x°.  However,  the  tabu  structure  keeps  the  escape  procedure  from 
converging  to  the  same  solution  each  time  and  we  do  not  continue  the  escape 
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procedure  for  very  long.  We  minimize  with  respect  to  Equation  6.1  until  a  local 
optima  is  reached  or  MAX_ESCAPE_ITER  iterations  have  been  performed. 

6.1.7  Cutting  Planes 

Since  we  are  using  xLP*  to  guide  the  search,  it  makes  sense  to  consider  adding 
cutting  planes  to  the  LP  to  move  jclp*  closer  to  the  optimal  IP  solution  and  hopefully 

improve  our  guide.  There  are  many  types  of  cutting  planes  each  with  their  own 
advantages  and  disadvantages  depending  on  the  problem  being  solved.  Instead  of 
developing  an  algorithm  to  generate  and  implement  of  these  cuts,  a  massive 
undertaking,  we  simply  let  CPLEX  apply  the  cuts. 

CPLEX  provides  several  callback  routines  that  allow  an  algorithm  to  interact 
with  CPLEX’ s  solvers  as  they  are  executing.  The  user  provides  CPLEX  the  address 
of  a  subroutine  to  be  called  when  a  specific  event  occurs.  When  the  user’s  routine  is 
completed,  control  is  passed  back  to  CPLEX  (unless  the  user  terminates  it). 

One  such  callback  for  the  mixed  integer  program  optimizer  (mipopt)  is  the 
heuristic  callback  which  calls  the  user’s  specified  routine  at  every  viable  node 
(feasible,  not  fathomed)  in  the  branch  and  bound  tree.  The  solution  to  the  current  LP 
relaxation,  xrLP* ,  as  well  as  pointers  that  can  be  used  to  access  the  other  information 

from  the  current  LP  is  passed  to  the  user’s  routine.  At  node  0  in  mipopt,  CPLEX 
adds  cuts  and  solves  the  node  LP  repeatedly  until  the  cuts  no  longer  make  significant 
progress.  CPLEX  invokes  the  heuristic  callback  after  solving  each  node  LP. 
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We  call  CPLEX’s  mipopt  to  solve  our  problem.  Cuts  added  at  node  0  are 
global  cuts  applying  to  the  whole  problem.  We  retrieve  the  new  solution,  basis  status, 
and  reduced  costs  after  every  node  0  LP.  Once  node  0  is  completed,  we  use  the  final 
node  LP  as  our  xLP*  and  terminate  CPLEX.  Which  JtLP*  we  used  affects  the  starting 
point,  the  escape  procedure,  and  bound  strengthening.  We  test  our  algorithm  using 
the  xLP*  with  and  without  cuts  to  determine  if  the  cuts  do  indeed  improve 
performance. 

6.1.8  Bound  Strengthening 

We  employ  the  bound  strengthening  in  our  algorithm  as  described  in  Section 
5.3.4.  Given  an  objective  function  value  zgmp,  the  bounds  on  Xj  are 

Ibj  =min^lbj,ub-(^zGMP/cj ^  +  for  upper  bound  non-basic  variable  j  and 

ubj  =  mi n  ( ubj ,  /b  -t-  ^ zCMP /cj  J )  for  lower  bound  non-basic  variable  j.  The  bounds  are 

updated  whenever  a  new  best  solution  is  discovered. 

CPLEX  also  performs  bound  strengthening  as  part  of  the  pre-processing  for 
each  node.  We  can  use  a  callback  function  as  described  in  the  previous  section  to 
retrieve  the  strengthened  bounds  from  node  0  in  the  mipopt  algorithm.  The  details  of 
the  bound  strengthening  performed  by  CPLEX  are  not  provided.  The  reduced  cost 
based  bound  strengthening  is  implemented  in  all  versions  of  GTTS-GMP  tested.  The 
algorithms  are  tested  with  and  without  the  CPLEX  bounds  to  determine  their  effect 
on  performance. 
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6.1.9  Strategic  Oscillation 

Strategic  oscillation  is  allowing  the  search  to  exit  the  feasible  region  at 
strategic  points  then  driving  it  back  inside  in  hopefully  a  new  area  of  the  feasible 
solution  space.  We  maintain  bounds  feasibility  of  the  decision  variables  throughout 
the  algorithm.  However,  our  initial  solution  may  be  infeasible  in  terms  of  the 
constraint  set.  We  attempt  to  reduce  infeasibility  at  each  iteration  until  the  feasible 
region  is  reached.  Once  feasible,  we  allow  only  feasible  moves.  If  the  neighborhood 
does  not  contain  a  feasible  non-tabu  move  the  best  infeasible  non-tabu  move  is 
chosen. 

We  incorporate  strategic  oscillation  in  both  escape  procedures.  In  the  random 
restart  strategy,  the  random  starting  solution  is  most  likely  infeasible.  We  then  move 
towards  the  feasible  region  as  described  above.  In  the  surrogate  objective  function 
strategy,  we  completely  ignore  feasibility  of  the  constraint  set  as  we  move  towards 
xLf'  allowing  the  search  to  exit  the  feasible  region  if  it  improves  the  objective. 

6.1.10  Stopping  Criteria 

Since  we  are  using  a  search  heuristic,  it  is  unlikely  we  will  be  able  to  prove 
that  a  global  optimal  solution  has  been  found.  Since  some  of  our  basic  variables  are 
fractional  we  know  at  least  one  of  the  non-basic  variables  must  change  by  at  least  1. 
Therefore,  z'LP  -min(cj)  is  an  upper  bound  on  a  global  optimal  solution  for  the  ILP. 
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If  we  find  a  feasible  solution  that  equals  this  bound  then  we  can  terminate  with  a 
global  optimal  solution.  Otherwise,  we  continue  to  search  until  a  maximum  time 
limit  or  number  of  iterations  has  been  reached. 


6.2  Computational  Results 


6.2.1  Test  Cases 

We  use  set  covering  and  multi-dimensional  knapsack  problem  sets  to  test  our 
algorithm.  These  problems  contain  a  diverse  set  of  characteristics.  SCPs  are 
minimization  problems  and  MDKPs  are  maximization  problems.  SCPs  are  binary 
while  MDKPs  can  be  binary  or  general  integer  programming  problems.  We  solve 
both  versions  here.  SCPs  have  relatively  sparse  constraint  matrices  consisting  of  all 
Is  and  Os  while  MDKPs  have  relatively  dense  matrices  containing  almost  any  values. 


problem  set 

number  of 
problems 

number  of 
columns 

number  of 

rows 

density 

4 

10 

1000 

%i§»»inSSji 

2% 

5 

2000 

2% 

6 

5 

1000 

200 

5% 

A 

5 

3000 

300 

2% 

B 

5 

3000 

300 

5% 

C 

5 

4000 

400 

2% 

D 

5 

4000 

400 

5% 

E 

5 

500 

50 

20% 

NRE 

5 

5000 

500 

10% 

NRF 

5 

5000 

500 

20% 

NRG 

5 

10000 

1000 

2% 

NRH 

5 

10000 

1000 

5% 

Table  6.1  Characteristics  for  SCP  sets 
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The  benchmark  problems  solved  were  obtained  from  Beasley’s  OR-Library 
(Beasley  1990).  For  set  covering,  problem  sets  4-6  originally  appeared  in  (Balas  and 


Ho  1980),  problem  sets  A-E  appeared  in  (Beasley  1987)  and  problem  sets  NRE-NRH 


appeared  in  (Beasley  1990).  All  problems  were  randomly  generated  based  on  the 
strategy  of  (Balas  and  Ho  1980).  All  of  the  problem  sets  except  E  were  generated  as 
weighted  SCPs.  Problem  set  E  is  a  unicost  SCP.  The  characteristics  of  the  SCP  sets 
are  shown  in  Table  6.1. 

For  knapsack,  problem  set  1  appeared  in  (Petersen  1967),  problem  sets  CB1- 
CB9  appeared  in  (Chu  and  Beasley  1998).  All  problems  were  generated  as  binary 


problems;  however,  we  solve  each  problem  as  binary  and  as  general  integer.  Problem 
set  1  contains  small  easy  to  solve  problems.  Problem  sets  CB1-CB9  get  progressively 
more  difficult.  There  are  30  problems  in  each  set.  We  solve  the  first  problem  in  sets 
CB1-CB8  and  the  first  10  problems  in  CB9.  The  characteristics  for  the  MDKP  sets 
are  shown  in  Table  6.2. 


problem  set 

number  of 
columns 

number  of 

rows 

binary 

1 

7 

6-50 

5-10 

n 

lb 

7 

5-10 

y 

cbl  l-cb81 

8 

100-500 

5-30 

n 

cbllb-cb81b 

8 

100-500 

5-30 

y 

cb91-cb910 

10 

500 

30 

n 

cb91b-cb910b 

10 

500 

30 

y 

Table  6.2  Characteristics  for  MDKP  sets 


6.2.2  Test  Procedures 

All  tests  were  performed  on  Dell  Precision  530  Workstations  running  SuSE 


Linux  with  two  1.8GHz  Pentium  Xeon  processors  and  1GB  of  RAM.  The  machines 
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are  multi-user  platforms.  An  attempt  was  made  to  find  machines  that  were  not  too 
busy,  but  as  some  problems  ran  for  at  least  an  hour,  CPU  usage  surely  fluctuated 
during  processing.  Each  problem  was  solved  using  CPLEX  9.0  and  GTTS-GMP  with 
a  time  limit  of  3600  seconds  for  either  method.  GTTS-GMP  was  coded  in  C. 

6.2.3  CPLEX  9.0 

To  provide  a  high  benchmark,  we  compare  our  results  to  CPLEX  version  9.0. 
CPLEX  uses  a  very  sophisticated  branch  and  cut  algorithm  to  solve  MILPs.  The 
algorithm  is  further  aided  by  two  heuristics.  The  first  attempts  to  create  a  feasible 
solution  from  the  fractional  solution  at  the  node.  The  second  attempts  to  improve  the 
incumbent  integer  solution  through  a  neighborhood  search  (ILOG  2003). 

It  would  be  naive  to  think  that  we  could  outperform  CPLEX  with  a  general 
ILP  algorithm.  CPLEX  was  developed  with  millions  of  dollars  and  man-centuries  of 
effort.  CPLEX  achieves  the  optimal  or  best  known  solutions  for  nearly  every  test 
problem.  We  simply  hope  to  perform  well  against  this  benchmark  to  demonstrate  the 
potential  of  our  approach. 

6.2.4  Results 

In  our  testing,  we  found  that  the  dynamic  neighborhoods  performed  better 
than  the  1-step  neighborhoods  on  the  MDKP  problem  sets.  However,  we  found  the 
opposite  to  be  true  for  the  SCP  problems  sets.  This  may  be  a  factor  of  neighborhood 
size.  The  largest  MDKP  test  instance  has  500  columns  while  the  smallest  SCP  test 
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instance  has  1000  columns.  So  for  the  SCP  sets  the  1-step  neighborhoods  are 
significantly  larger  than  for  the  MDKP  sets. 

We  used  the  dynamic  neighborhood  option  for  the  MDKP  sets  and  the  1-step 
neighborhood  for  the  SCP  sets  for  all  algorithm  configurations.  The  tabu  tenure  also 
differs  between  problem  sets.  The  default  move  tenure  is  a  function  of  problem  size, 
n*0.1.  The  default  solution  tabu  tenure  is  smaller  when  dynamic  neighborhoods  are 
used.  It  makes  sense  to  revisit  solutions  more  frequently  when  the  neighborhood  has 
changed.  Since  we  are  using  reactive  tabu  search,  both  tabu  tenures  change 
throughout  the  search. 

We  next  examine  the  surrogate  objective  function  escape  strategy.  The  results 
for  the  SCP  sets  are  presented  in  Table  6.3.  The  results  for  the  MDKP  problem  sets 
are  presented  in  Table  6.4.  We  tested  each  problem  with  cuts  and  bounds  from 
CPLEX,  with  cuts  only,  with  bounds  only,  and  without  cuts  or  bounds. 

With  the  exception  of  the  smaller  problems  adding  cuts  and/or  bounds  from 
CPLEX  improves  the  performance  on  the  SCP  sets.  The  bounds  from  CPLEX  for 
SCPs  are  generally  stronger  than  those  derived  using  reduced  costs.  Also  the 
algorithm  has  these  bounds  at  the  start  and  does  not  need  to  wait  for  new  solutions  to 
generate  them.  The  bounds  seem  to  affect  performance  more  than  the  cuts.  All 
things  considered,  GTTS-GMP  with  cuts  and  bounds  performs  best. 

For  the  MDKP  sets  the  effect  of  bounds  and  cuts  is  less  significant.  For  the 
general  integer  problems,  the  bounds  from  reduced  cuts  quickly  dominate  those  from 
CPLEX.  In  either  case  the  bounds  from  CPLEX  are  not  as  strong  as  they  are  for  the 
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SCP  sets.  In  general,  the  cuts  do  not  have  as  much  of  an  impact  on  the  algorithm  as 
we  anticipated.  However,  we  still  feel  that  GTTS-GMP  with  cuts  and  bounds 
performs  best. 


problem 

set 


optimal  or 
*best  known 
average 


CPLEX 
average 
(time  sec) 


GTTS-GMP  average  -  (time  in  seconds)  -  %  gap 


cuts  & 
bounds 


cuts  only 


bounds 

only 


none 


510 


510 

(0) 


511.2 

(10.2) 

0.22% 


512.3 

(3.7) 

0.45% 


511.3 

(0.5) 

0.25% 


510.4 

(0.9) 

0.08% 


257.2 


257.2 

(0.3) 


259.2 

(8) 

0.78% 


259 

(12) 

0.70% 


259.2 

(6.9) 

0.78% 


258.8 

(6.5) 

0.62% 


144.2 


144.2 

(0) 


147.8 

(41.6) 

2.5% 


148.2 
(  32.2 ) 
2.77% 


146.6 

(6.8) 

1.66% 


148.2 

(5.4) 

2.77% 


241.4 


241.4 

(1.2) 


245.6 

( 54.2 ) 
1.74% 


246.8 

(69) 

2.24% 


246.4 

( 73.2 ) 
2.07% 


246.6 

(19.2) 

2.15% 


B 


75.2 


75.2 

(1.8) 


76.8 
( 56.8  ) 
2.13% 


77.8 
( 44.8  ) 
3.46% 


77 

(61.2) 

2.39% 


77.4 

(29.4) 

2.93% 


224.6 


224.6 

(1.8) 


233 

(88) 

3.74% 


233.4 
(  107.2 ) 
3.92% 


232.4 

(41.4) 

3.47% 


232.6 

(131) 

3.56% 


D 


64.2 


64.2 

(8) 


66.8 

(65.8) 

4.05% 


66.4 
(  117.2) 
3.43% 


67.4 

(55) 

4.98% 


67 

(138) 

4.36% 


5 

(0) 


5 

(0) 

0.0% 


5 

(0) 

0.0% 


5 

(0) 

0.0% 


5 

(0) 

0.0% 


NRE 


28.4 


28.4 
(  390.2 ) 


29.4 

(147) 

3.52% 


29.6 

(291.4) 

4.23% 


29 

(  785.8  ) 
2.11% 


29.8 
( 56.4 ) 
4.93% 


NRF 


14 


14 

(721) 


14.2 

(14) 

1.43% 


14.8 

(  105.2) 
5.71% 


14 

(621.2) 

0.0% 


15 

(8.8) 

7.14% 


NRG 


*166.4 


166.8 
( 625.6 ) 


180 

(479) 

8.17% 


181.6 

(1180) 

9.13% 


180 

( 430.6 ) 
8.17% 


180.6 
(  1522 ) 
8.53% 


NRH 


*59.6 


61.2 

(  1361.2) 


65.8 

(183) 

10.40% 


66 

( 798.6 ) 
10.74% 


66 

( 226.8 ) 
10.74% 


68 

(208) 

14.09% 


Table  6.3  SCP  results  for  surrogate  objective  function  escape  strategy 


problem 

set 

optimal  or 
*best  known 
average 

CPLEX 
average 
(time  sec) 

GTTS-GMP  avera 

ge  (time  in  seconds) 

cuts  & 
bounds 

cuts  only 

none 

1 

13604.41 

13604.41 

(0) 

13604.27 

(0) 

0.0% 

13603.56 

(0.3) 

0.01% 

13603.56 

(0.6) 

0.01% 

13604.27 

(0.1) 

0.0% 

lb 

8885.16 

8885.16 

(0) 

8883.16 

(0.1) 

0.02% 

8883.16 

(1.3) 

0.02% 

8885.16 

(0) 

0.0% 

8885.16 

(1.4) 

0.0% 

cbll-cbl8 

*78337.5 

78337.5 
(  166.5 ) 

78173.63 

(119.5) 

0.21% 

78157.63 

(311.3) 

0.23% 

78179.75 
( 430.3 ) 
0.20% 

78179.63 

(295) 

0.20% 

cbllb- 

cbl8b 

*60330.13 

60330.13 
( 441.38 ) 

60192.5 

(515.9) 

0.23% 

60202 
( 849.5 ) 
0.21% 

cb91-cb910 

*125682.8 

125682.8 
(  1773.5 ) 

125194.6 
( 1118.6) 
0.39% 

125165.4 
( 1391.5) 
0.41% 

125163.2 
(  1020.9) 
0.41% 

125154.1 

(1211) 

0.42% 

cb91b- 

cb910b 

*115553.9 

115553.9 
(  1736.9) 

115229.4 
(  1348.6) 
0.28% 

115227.5 

(918.2) 

0.28% 

115232.8 
(  1122.5) 
0.28% 

Table  6.4  MDKP  results  for  surrogate  objective  function  escape  strategy 


Next  we  examine  the  random  restart  escape  strategy.  The  results  for  the  SCP 
sets  are  presented  in  Tables  6.5.  The  results  for  the  MDKP  problem  sets  are 
presented  in  Tables  6.6.  We  tested  each  problem  with  cuts  and  bounds  from  CPLEX, 
with  cuts  only,  with  bounds  only,  and  without  cuts  or  bounds.  The  results  are  the 
average  of  three  runs  for  each  configuration. 

With  the  exception  of  set  1  adding  cuts  and/or  bounds  from  CPLEX  improves 
the  performance  on  the  SCP  sets.  The  bounds  have  a  much  greater  impact  on 
performance  than  the  cuts.  One  possible  reason  for  this  is  that  the  cuts  reduce  the 
fractional  parts  of  xLP"  which  then  limits  the  variability  of  the  random  restart 

solutions.  On  average  we  still  feel  that  GTTS-GMP  with  cuts  and  bounds  performs 
best. 
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Again  for  the  MDKP  sets  the  effect  of  bounds  and  cuts  is  less  significant  for 
the  same  reasons  as  above.  GTTS-GMP  with  cuts  and  bounds  also  performs  best  on 
he  MDKP  sets. 


problem 

set 


optimal  or 
*best  known 
average 


CPLEX 
average 
(time  sec) 


GTTS-GMP  average  -  (time  in  seconds)  -  %  gap 


cuts  & 
bounds 


cuts  only 


bounds 

only 


none 


510 


510 

(0) 


512.37 

(6.3) 

0.46% 


511.9 
(  10.1  ) 
0.37% 


514.06 
( 3.23 ) 
0.80% 


510.2 
( 12.5 ) 
0.04% 


257.2 


257.2 

(0.3) 


258 
(  17.5  ) 
0.31% 


258.33 

(18) 

0.44% 


258.6 
( 24.4 ) 
0.54% 


258.17 

(15.9) 

0.38% 


144.2 


144.2 

(0) 


145.13 
(  19.4 ) 
0.65% 


145.6 

(41.9) 

0.97% 


145 

(31.5) 

0.55% 


145.53 
( 36.7 ) 
0.92% 


241.4 


241.4 

(1.2) 


243.6 

( 79.9 ) 
0.91% 


243.93 
(  103.3) 
1.05% 


244.2 

(84.1) 

1.16% 


244.8 
( 101.4) 
1.41% 


B 


75.2 


75.2 

(1-8) 


75.87 

( 57.9 ) 
0.89% 


76.2 
( 60.1  ) 
1.33% 


75.8 
( 105.2) 
0.80% 


76.27 
( 85.5  ) 
1.42% 


224.6 


224.6 

(1.8) 


230.53 

(108) 

2.64% 


230.53 
( 133.5) 
2.64% 


230.33 

(78.9) 

2.55% 


231.33 
( 124.1 ) 
3.00% 


D 


64.2 


64.2 

(8) 


65.13 
( 92.3  ) 
1.45% 


65.73 

(111.7) 

2.39% 


65.47 
( 88.5 ) 
1.97% 


65.87 

(80.1) 

2.60% 


5 

(0) 


5 

(0) 

0.0% 


5 

(0) 

0.0% 


5 

(0) 

0.0% 


5 

(0) 

0.0% 


NRE 


28.4 


28.4 
( 390.2  ) 


28.6 
( 433.3  ) 
0.70% 


29 

(  146.67 ) 
2.11% 


28.53 

(930) 

0.47% 


28.73 
( 654.5  ) 
1.17% 


NRF 


14 


14 

(721) 


14.07 

(341.1) 

0.48% 


14.2 
( 62.8  ) 
1.43% 


14.13 
(239.1  ) 
0.95% 


14.07 
(345.1  ) 
0.48% 


NRG 


*166.4 


166.8 
( 625.6 ) 


179.07 
( 605.3  ) 
7.61% 


178.93 
(612.1  ) 
7.53% 


178.5 

(1132.3) 

7.29% 


178.4 

(391.3) 

7.21% 


NRH 


*59.6 


61.2 

( 1361.2) 


63.67 
(  1078.1  ) 
6.82% 


63.33 
(  1015.8) 
6.26% 


63.73 

( 564.3 ) 
6.94% 


63.73 
(  1241.9) 
6.94% 


Table  6.5  SCP  results  for  random  restart  escape  strategy 
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problem 

set 

optimal  or 
*best  known 
average 

CPLEX 
average 
(time  sec) 

GTTS-GMP  avera 

ge  (time  in  seconds) 

cuts  & 
bounds 

cuts  only 

none 

1 

13604.41 

13604.41 

(0) 

13604.03 

(0.1) 

0.0% 

13604.03 

(0.3) 

0.0% 

13603.27 

(0.4) 

0.00% 

13600.6 

(0.8) 

0.03% 

lb 

8885.16 

8885.16 

(0) 

8884.49 

(0.2) 

0.01% 

8883.16 

(0.1) 

0.02% 

8883.59 

(1.1) 

0.02% 

8877.97 

(0.2) 

0.08% 

cbll-cbl8 

*78337.5 

78337.5 
(  166.5  ) 

78199.42 
(  222.6 ) 
0.18% 

78205.83 
( 307.8  ) 
0.17% 

78198.92 

(341.1) 

0.18% 

78200.96 
( 540.9 ) 
0.17% 

cbllb- 

cbl8b 

*60330.13 

60330.13 

(441.38) 

60225.5 
( 555.6 ) 
0.17% 

60229.46 

(513.6) 

0.17% 

60205.5 
( 452.9 ) 
0.21% 

60198.21 
( 669.83  ) 
0.22% 

cb91-cb910 

*125682.8 

125682.8 
(  1773.5  ) 

125379.7 
(  1685.8 ) 
0.24% 

125382.7 
(  1374.7 ) 
0.24% 

125398.3 
(  1384.9) 
0.23% 

125372.8 
(  1782.5  ) 
0.25% 

cb91b- 

cb910b 

*115553.9 

115553.9 
(  1736.9 ) 

115361.9 
( 1615.9) 
0.17% 

115323.8 
(  1311.4) 
0.20% 

115336.6 
(  1511.5) 
0.19% 

115313.1 
( 1781.1 ) 
0.21% 

Table  6.6  MDKP  results  for  random  restart  escape  strategy 


The  random  restart  escape  strategy  is  clearly  superior  to  the  surrogate 
objective  function.  Although  both  strategies  are  able  to  escape  from  cycling,  the 
surrogate  objective  function  drives  the  search  towards  the  same  point,  x0,  every  time. 
While  the  tabu  structure  stops  the  procedure  from  visiting  the  same  solution  from  one 
escape  to  the  next,  it  still  does  not  provide  the  diversification  that  we  get  from  the 
random  restart  strategy. 


6.3  Super  Optimal  Solutions 

Our  algorithm  oscillates  in  and  out  of  the  feasible  region  during  the  search. 


Often  high  quality  solutions  can  be  found  with  very  minor  violations  in  feasibility. 


We  call  such  solutions  super-optimal  (Carlton  and  Barnes  1996).  The  constraints  in 
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our  models  are  often  the  result  of  some  simplifying  assumptions  and  are  rarely 


pristine.  The  decision  maker  may  want  to  evaluate  these  super-optimal  solutions  to 
determine  if  they  could  accept  the  violations  given  the  value  of  the  solution.  Tables 
6.7  and  6.8  contain  examples  of  the  super-optimal  solutions  found.  For  MDKP  the 
violation  is  normalized  by  the  right  hand  side.  For  example,  if  the  constraint  is  <  500 


and  the  value  is  525  then  the  violation  is  25/500  or  0.05. 


problem 

optimal 

super-optimal 

violation 

44 

497 

470 

1 

49 

641 

630 

1 

58 

288 

265 

1 

64 

131 

123 

1 

Table  6.7  Some  super-optimal  solutions  from  SCP  sets 


problem 

optimal 

super-optimal 

violation 

12 

10970.9 

1 1082.7 

0.004 

12b 

8706.1 

8886.2 

0.009 

cb21 

93127 

93216 

0.0005 

cb31 

175856 

176082 

0.0009 

cb71b 

21946 

22008 

0.007 

Table  6.8  Some  super-optimal  solutions  from  MDKP  sets 


6.4  Conclusion 

Our  GTTS-GMP  algorithm  performs  well,  finding  solutions  well  within  5%  of 
the  best  known  for  all  but  2  problem  sets.  Many  of  these  techniques,  rounding  the  LP 
relaxation,  bounding  by  reduced  costs,  using  cuts,  etc.,  are  universal  and  can  be 
applied  to  tailored  algorithms  designed  for  specific  problems. 

For  our  test,  the  random  restart  strategy  with  cuts  and  variable  bounds  from 
CPLEX  performed  the  best  overall.  More  research  needs  to  be  done  to  improve  the 
performance  of  the  general  algorithm.  However,  even  if  a  general  algorithm  is  used  it 


100 


is  clear  that  tuning  parameters,  such  as  dynamic  neighborhood  and  tenure,  for  a 
specific  problem  type  will  improve  performance. 
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Chapter  7  -  Conclusion  &  Future  Research 


7.1  Conclusions 

In  this  work  we  explored  the  use  of  group  theory  in  metaheuristic  search 
methods  for  partitioning  integer  linear  programs  (ILPs).  To  the  best  of  our 
knowledge,  this  is  the  first  time  local  search  methods  for  these  types  of  problems 
have  been  viewed  through  this  framework.  Past  efforts  have  demonstrated  the  power 
of  group  theory  in  metaheuristic  search  methods  for  problems  with  an  ordering 
component. 

7.1.1  Partitioning  into  Orbits 

Using  group  theory  we  defined  procedures  for  partitioning  the  solution  space 
into  orbits.  By  clustering  the  variables,  we  are  able  to  create  “good”  and  “bad” 
orbits.  We  developed  neighborhoods  to  explore  the  individual  orbits  and  transition 
between  them.  We  also  developed  methods  for  bounding  the  solutions  in  each  orbit 
so  that  the  may  be  discarded  as  infeasible  or  dominated  by  the  incumbent  solution. 

We  tested  these  techniques  by  developing  a  group  theoretic  tabu  search 
algorithm  to  solve  the  unicost  set  covering  problem,  the  GTTS-USCP.  Our  variable 
clustering  was  based  on  the  reduced  costs  of  the  LP  relaxation  of  the  problem.  The 
use  of  variable  clustering  and  group  theory  allowed  our  algorithm  to  intensify  the 
search  in  the  areas  of  the  solution  space  believed  to  contain  good  solutions.  The 
orbits  kept  the  search  contained  in  these  areas  and  the  clusters  worked  as  an  enhanced 
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candidate  list,  reducing  the  total  number  of  moves  in  the  neighborhood  while  still 
retaining  the  “good”  moves.  These  techniques  proved  very  effective  for  the  unicost 
SCP  discovering  46  new  best  known  solutions  for  the  set  of  65  benchmark  problems 
and  outperformed  CPLEX  in  both  solution  quality  and  speed. 

7.1.2  The  Group  Minimization  Problem 

Next  we  looked  at  the  general  ILP  by  examining  its  associated  group 
minimization  problem  (GMP).  We  demonstrated  why  valid  column  reduction  for  the 
GMP  is  not  valid  in  terms  of  the  original  ILP.  We  examined  the  comer  polyhedron  in 
the  space  of  the  set  of  non-basic  variables.  We  proved  that  the  density  of  the  integer 
points  in  non-basic  variable  space  that  yield  an  all-integer  basis  is  bounded  by  the 
determinant  of  the  basis.  We  defined  new  search  neighborhoods,  identity  move 
neighborhoods,  to  traverse  the  sub-lattice  formed  by  these  points.  Finally,  we 
developed  procedures  for  strengthening  the  bounds  on  the  non-basic  variables  using 
the  reduced  cost  from  the  optimal  solution  to  the  LP  relaxation.  These  bounds  reduce 
the  size  of  the  comer  polyhedron  and  the  solution  space  of  the  GMP. 

Based  on  these  results  we  developed  a  GTTS  algorithm  for  the  GMP,  GTTS- 
GMP.  Since  a  GMP  can  be  formed  for  any  ILP,  our  algorithm  solves  the  general  ILP. 
We  are  able  to  add  cuts  and  additional  variable  bounds  from  CPLEX  at  the  beginning 
of  our  algorithm.  The  algorithm  performs  well  against  a  diverse  set  of  test  problems. 
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7.2  Future  Research 


7.2.1  Other  Clustering  Techniques 

For  the  unicost  SCP  problem,  the  LP  solution  provided  a  good  profile  of 
quality  ILP  solutions  and  provided  a  good  basis  for  variable  clustering.  This  is  not 
likely  to  be  true  for  all  problems  types.  In  addition,  for  problems  where  solving  the 
LP  relaxation  is  not  reasonable  other  methods  of  clustering  will  need  to  be  developed. 
Heuristic  methods,  such  as  benefit/cost  ratio  for  knapsack  problems,  may  provide  a 
more  efficient  and  accurate  clustering  and  should  be  explored. 

Another  possible  future  enhancement  to  the  partitioning  algorithm  is  the  use 
of  composite  moves.  By  using  only  single  transpositions,  we  only  implemented 
simple  swap  moves  in  our  move  neighborhoods.  It  might  be  worthwhile  to  create 
composite  moves  by  combining  transpositions  from  each  cluster. 

7.2.2  Embed  GTTS-GMP  in  Branch  &  Cut 

GTTS-GMP  could  be  embedded  within  a  branch  &  cut  program  such  as 
CPLEX.  The  branch  &  cut  algorithm  could  execute  a  short  version  of  the  GTTS- 
GMP  algorithm  at  each  viable  node  to  attempt  to  find  a  new  feasible  incumbent 
solution.  The  local  bounds  at  the  node  provide  diversification  for  the  GTTS-GMP 
algorithm  and  the  solutions  found  by  GTTS-GMP  help  the  branch  &  cut  algorithm  to 
fathom  more  nodes.  We  feel  this  concept  holds  a  lot  of  promise  and  should  be 
explored  further. 
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7.2.3  Equality  Constraints 

The  method  for  creating  identity  moves  presented  in  Section  5.3.2  does  not 
consider  equality  constraints.  One  approach  may  be  to  replace  each  equality 
constraint  with  a  >  and  <  constraint  and  proceed  as  normal.  This  may  present 
difficultly  in  maintaining  a  feasible  solution  during  the  search.  Another  approach 
may  be  to  create  the  1-step  neighborhood  ignoring  the  equality  constraints  and 
penalize  the  search  based  on  the  violation  of  these  constraints.  It  may  also  be 
possible  to  use  the  coefficients  in  the  equality  constraints  to  combine  1-step  moves 
into  composite  moves  that  maintain  the  feasibility  of  these  constraints.  An  efficient 
and  effective  method  for  dealing  with  these  constraints  is  another  area  of  research. 

7.2.4  Mixed  Integer  Linear  Programs 

The  GTTS-GMP  is  for  all  integer  problems.  How  can  we  apply  these 
techniques  to  mixed  integer  problems?  There  is  a  mixed  integer  version  of  the  GMP 
developed  by  Araoz  (1973).  The  columns  of  his  problem  form  an  abelian  semigroup. 
A  semigroup  lacks  the  inverse  and  identity  properties  of  a  group.  Perhaps  the 
semigroup  minimization  problem  can  help  us  develop  techniques  for  the  mixed 
integer  case. 

The  research  documented  here  was  supported  by  a  grant  for  the  Air  Force  Office  of 
Scientific  research. 
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